Index: /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c	(revision 1593)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c	(revision 1594)
@@ -47,81 +47,5 @@
 }
 
-int find_qmin_and_qmax(double dq0, double dq1, double dq2, double *qmin, double *qmax){
-  //Considering the centroid of an FV triangle and the vertices of its auxiliary triangle, find
-  //qmin=min(q)-qc and qmax=max(q)-qc, where min(q) and max(q) are respectively min and max over the
-  //four values (at the centroid of the FV triangle and the auxiliary triangle vertices),
-  //and qc is the centroid
-  //dq0=q(vertex0)-q(centroid of FV triangle)
-  //dq1=q(vertex1)-q(vertex0)
-  //dq2=q(vertex2)-q(vertex0)
-  if (dq0>=0.0){
-    if (dq1>=dq2){
-      if (dq1>=0.0)
-	*qmax=dq0+dq1;
-      else
-	*qmax=dq0;
-      if ((*qmin=dq0+dq2)<0)
-	;//qmin is already set to correct value
-      else
-	*qmin=0.0;
-    }
-    else{//dq1<dq2
-      if (dq2>0)
-	*qmax=dq0+dq2;
-      else
-	*qmax=dq0;
-      if ((*qmin=dq0+dq1)<0)
-	;//qmin is the correct value
-      else
-	*qmin=0.0;
-    }
-  }
-  else{//dq0<0
-    if (dq1<=dq2){
-      if (dq1<0.0)
-	*qmin=dq0+dq1;
-      else
-	*qmin=dq0;
-      if ((*qmax=dq0+dq2)>0.0)
-	;//qmax is already set to the correct value
-      else
-	*qmax=0.0;
-    }
-    else{//dq1>dq2
-      if (dq2<0.0)
-	*qmin=dq0+dq2;
-      else
-	*qmin=dq0;
-      if ((*qmax=dq0+dq1)>0.0)
-	;//qmax is already set to the correct value
-      else
-	*qmax=0.0;
-    }
-  }
-  return 0;
-}
-
-int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
-  //given provisional jumps dqv from the FV triangle centroid to its vertices and
-  //jumps qmin (qmax) between the centroid of the FV triangle and the
-  //minimum (maximum) of the values at the centroid of the FV triangle and the auxiliary triangle vertices,
-  //calculate a multiplicative factor phi by which the provisional vertex jumps are to be limited
-  int i;
-  double r=1000.0, r0=1.0, phi=1.0;
-  static double TINY = 1.0e-100;//to avoid machine accuracy problems.
-  //Any provisional jump with magnitude < TINY does not contribute to the limiting process.
-  for (i=0;i<3;i++){
-    if (dqv[i]<-TINY)
-      r0=qmin/dqv[i];
-    if (dqv[i]>TINY)
-      r0=qmax/dqv[i];
-    r=min(r0,r);
-    //
-  }
-  phi=min(r*beta_w,1.0);
-  for (i=0;i<3;i++)
-    dqv[i]=dqv[i]*phi;
-  return 0;
-}
+
 
 // Computational function for flux computation (using stage w=z+h)
@@ -246,6 +170,5 @@
       if (h >= eps) {
 	S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
-	//S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
-	S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
+	S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
 	//S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
 
@@ -503,292 +426,4 @@
 }
 
-PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
-  /*Compute the vertex values based on a linear reconstruction on each triangle
-    These values are calculated as follows:
-    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
-    formed by the centroids of its three neighbours.
-    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
-    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
-    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
-    original triangle has gradient (a,b).
-    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
-    find_qmin_and_qmax and limit_gradient
-
-    Python call:
-    extrapolate_second_order_sw(domain.surrogate_neighbours,
-                                domain.number_of_boundaries
-                                domain.centroid_coordinates,
-                                Stage.centroid_values
-                                Xmom.centroid_values
-                                Ymom.centroid_values
-                                domain.vertex_coordinates,
-                                Stage.vertex_values,
-                                Xmom.vertex_values,
-                                Ymom.vertex_values)
-
-    Post conditions:
-            The vertices of each triangle have values from a limited linear reconstruction
-	    based on centroid values
-
-  */
-  PyArrayObject *surrogate_neighbours,
-    *number_of_boundaries,
-    *centroid_coordinates,
-    *stage_centroid_values,
-    *xmom_centroid_values,
-    *ymom_centroid_values,
-    *vertex_coordinates,
-    *stage_vertex_values,
-    *xmom_vertex_values,
-    *ymom_vertex_values;
-  PyObject *domain, *Tmp;
-  //Local variables
-  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
-  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
-  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
-  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
-  double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
-  //by which these jumps are limited
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
-			&domain,
-			&surrogate_neighbours,
-			&number_of_boundaries,
-			&centroid_coordinates,
-			&stage_centroid_values,
-			&xmom_centroid_values,
-			&ymom_centroid_values,
-			&vertex_coordinates,
-			&stage_vertex_values,
-			&xmom_vertex_values,
-			&ymom_vertex_values)) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
-    return NULL;
-  }
-
-  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
-  Tmp = PyObject_GetAttrString(domain, "beta_w");
-  if (!Tmp)
-    return NULL;
-  beta_w = PyFloat_AsDouble(Tmp);
-  Py_DECREF(Tmp);
-  number_of_elements = stage_centroid_values -> dimensions[0];
-  for (k=0; k<number_of_elements; k++) {
-    k3=k*3;
-    k6=k*6;
-
-    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
-      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
-      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
-      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
-      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
-      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
-      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
-      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
-      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
-      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
-      continue;
-    }
-    else {//we will need centroid coordinates and vertex coordinates of the triangle
-      //get the vertex coordinates of the FV triangle
-      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
-      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
-      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
-      //get the centroid coordinates of the FV triangle
-      coord_index=2*k;
-      x=((double *)centroid_coordinates->data)[coord_index];
-      y=((double *)centroid_coordinates->data)[coord_index+1];
-      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
-      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
-      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
-    }
-    if (((long *)number_of_boundaries->data)[k]<=1){
-      //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
-      //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
-      k0=((long *)surrogate_neighbours->data)[k3];
-      k1=((long *)surrogate_neighbours->data)[k3+1];
-      k2=((long *)surrogate_neighbours->data)[k3+2];
-      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
-      coord_index=2*k0;
-      x0=((double *)centroid_coordinates->data)[coord_index];
-      y0=((double *)centroid_coordinates->data)[coord_index+1];
-      coord_index=2*k1;
-      x1=((double *)centroid_coordinates->data)[coord_index];
-      y1=((double *)centroid_coordinates->data)[coord_index+1];
-      coord_index=2*k2;
-      x2=((double *)centroid_coordinates->data)[coord_index];
-      y2=((double *)centroid_coordinates->data)[coord_index+1];
-      //store x- and y- differentials for the vertices of the auxiliary triangle
-      dx1=x1-x0; dx2=x2-x0;
-      dy1=y1-y0; dy2=y2-y0;
-      //calculate 2*area of the auxiliary triangle
-      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
-      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
-      if (area2<=0)
-	return NULL;
-
-      //### stage ###
-      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
-      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
-      //calculate differentials between the vertices of the auxiliary triangle
-      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
-      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
-      //calculate the gradient of stage on the auxiliary triangle
-      a = dy2*dq1 - dy1*dq2;
-      a /= area2;
-      b = dx1*dq2 - dx2*dq1;
-      b /= area2;
-      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
-      dqv[0]=a*dxv0+b*dyv0;
-      dqv[1]=a*dxv1+b*dyv1;
-      dqv[2]=a*dxv2+b*dyv2;
-      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
-      //and compute jumps from the centroid to the min and max
-      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
-
-      //### xmom ###
-      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
-      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
-      //calculate differentials between the vertices of the auxiliary triangle
-      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
-      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
-      //calculate the gradient of xmom on the auxiliary triangle
-      a = dy2*dq1 - dy1*dq2;
-      a /= area2;
-      b = dx1*dq2 - dx2*dq1;
-      b /= area2;
-      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
-      dqv[0]=a*dxv0+b*dyv0;
-      dqv[1]=a*dxv1+b*dyv1;
-      dqv[2]=a*dxv2+b*dyv2;
-      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
-      //and compute jumps from the centroid to the min and max
-      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
-
-      //### ymom ###
-      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
-      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
-      //calculate differentials between the vertices of the auxiliary triangle
-      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
-      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
-      //calculate the gradient of xmom on the auxiliary triangle
-      a = dy2*dq1 - dy1*dq2;
-      a /= area2;
-      b = dx1*dq2 - dx2*dq1;
-      b /= area2;
-      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
-      dqv[0]=a*dxv0+b*dyv0;
-      dqv[1]=a*dxv1+b*dyv1;
-      dqv[2]=a*dxv2+b*dyv2;
-      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
-      //and compute jumps from the centroid to the min and max
-      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
-    }//if (number_of_boundaries[k]<=1)
-    else{//number_of_boundaries==2
-      //one internal neighbour and gradient is in direction of the neighbour's centroid
-      //find the only internal neighbour
-      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
-	if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
-	  break;
-      }
-      if ((k2==k3+3))//if we didn't find an internal neighbour
-	return NULL;//error
-      k1=((long *)surrogate_neighbours->data)[k2];
-      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
-      coord_index=2*k1;
-      x1=((double *)centroid_coordinates->data)[coord_index];
-      y1=((double *)centroid_coordinates->data)[coord_index+1];
-      //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
-      dx1=x1-x; dy1=y1-y;
-      //set area2 as the square of the distance
-      area2=dx1*dx1+dy1*dy1;
-      //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
-      //respectively correspond to the x- and y- gradients of the conserved quantities
-      dx2=1.0/area2;
-      dy2=dx2*dy1;
-      dx2*=dx1;
-
-      //## stage ###
-      //compute differentials
-      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
-      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
-      a=dq1*dx2;
-      b=dq1*dy2;
-      //calculate provisional vertex jumps, to be limited
-      dqv[0]=a*dxv0+b*dyv0;
-      dqv[1]=a*dxv1+b*dyv1;
-      dqv[2]=a*dxv2+b*dyv2;
-      //now limit the jumps
-      if (dq1>=0.0){
-	qmin=0.0;
-	qmax=dq1;
-      }
-      else{
-	qmin=dq1;
-	qmax=0.0;
-      }
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
-
-      //## xmom ###
-      //compute differentials
-      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
-      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
-      a=dq1*dx2;
-      b=dq1*dy2;
-      //calculate provisional vertex jumps, to be limited
-      dqv[0]=a*dxv0+b*dyv0;
-      dqv[1]=a*dxv1+b*dyv1;
-      dqv[2]=a*dxv2+b*dyv2;
-      //now limit the jumps
-      if (dq1>=0.0){
-	qmin=0.0;
-	qmax=dq1;
-      }
-      else{
-	qmin=dq1;
-	qmax=0.0;
-      }
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
-
-      //## ymom ###
-      //compute differentials
-      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
-      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
-      a=dq1*dx2;
-      b=dq1*dy2;
-      //calculate provisional vertex jumps, to be limited
-      dqv[0]=a*dxv0+b*dyv0;
-      dqv[1]=a*dxv1+b*dyv1;
-      dqv[2]=a*dxv2+b*dyv2;
-      //now limit the jumps
-      if (dq1>=0.0){
-	qmin=0.0;
-	qmax=dq1;
-      }
-      else{
-	qmin=dq1;
-	qmax=0.0;
-      }
-      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
-      for (i=0;i<3;i++)
-	((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
-    }//else [number_of_boudaries==2]
-  }//for k=0 to number_of_elements-1
-  return Py_BuildValue("");
-}//extrapolate_second-order_sw
-
 PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
   //
@@ -848,4 +483,5 @@
 }
 
+
 PyObject *compute_fluxes(PyObject *self, PyObject *args) {
   /*Compute all fluxes and the timestep suitable for all volumes
@@ -866,5 +502,5 @@
     domain.timestep = compute_fluxes(timestep,
                                      domain.epsilon,
-                                     domain.g,
+				     domain.g,
                                      domain.neighbours,
                                      domain.neighbour_edges,
@@ -882,6 +518,5 @@
                                      Stage.explicit_update,
                                      Xmom.explicit_update,
-                                     Ymom.explicit_update,
-                                     already_computed_flux)
+                                     Ymom.explicit_update)
 
 
@@ -905,6 +540,5 @@
     *stage_explicit_update,
     *xmom_explicit_update,
-    *ymom_explicit_update,
-    *already_computed_flux;//tracks whether the flux across an edge has already been computed
+    *ymom_explicit_update;
 
 
@@ -912,13 +546,12 @@
   double timestep, max_speed, epsilon, g;
   double normal[2], ql[3], qr[3], zl, zr;
-  double edgeflux[3]; //Work arrays for summing up fluxes
-
-  int number_of_elements, k, i, m, n;
-  int ki, nm=0, ki2; //Index shorthands
-  static long call=1;
+  double flux[3], edgeflux[3]; //Work arrays for summing up fluxes
+
+  int number_of_elements, k, i, j, m, n;
+  int ki, nm, ki2; //Index shorthands
 
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOO",
+  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO",
 			&timestep,
 			&epsilon,
@@ -937,24 +570,20 @@
 			&stage_explicit_update,
 			&xmom_explicit_update,
-			&ymom_explicit_update,
-			&already_computed_flux)) {
+			&ymom_explicit_update)) {
     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     return NULL;
   }
+
   number_of_elements = stage_edge_values -> dimensions[0];
-  call++;//a static local variable to which already_computed_flux is compared
-  //set explicit_update to zero for all conserved_quantities.
-  //This assumes compute_fluxes called before forcing terms
+
+
   for (k=0; k<number_of_elements; k++) {
-    ((double *) stage_explicit_update -> data)[k]=0.0;
-    ((double *) xmom_explicit_update -> data)[k]=0.0;
-    ((double *) ymom_explicit_update -> data)[k]=0.0;
-  }
-  //Loop through neighbours and compute edge flux for each
-  for (k=0; k<number_of_elements; k++) {
+
+    //Reset work array
+    for (j=0; j<3; j++) flux[j] = 0.0;
+
+    //Loop through neighbours and compute edge flux for each
     for (i=0; i<3; i++) {
       ki = k*3+i;
-      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
-	continue;
       ql[0] = ((double *) stage_edge_values -> data)[ki];
       ql[1] = ((double *) xmom_edge_values -> data)[ki];
@@ -972,4 +601,5 @@
       } else {
 	m = ((long *) neighbour_edges -> data)[ki];
+
 	nm = n*3+m;
 	qr[0] = ((double *) stage_edge_values -> data)[nm];
@@ -978,4 +608,9 @@
 	zr =    ((double *) bed_edge_values -> data)[nm];
       }
+
+      //printf("%d %d [%d] (%d, %d): %.2f %.2f %.2f\n", k, i, ki, n, m,
+      //	     ql[0], ql[1], ql[2]);
+
+
       // Outward pointing normal vector
       // normal = domain.normals[k, 2*i:2*i+2]
@@ -983,4 +618,5 @@
       normal[0] = ((double *) normals -> data)[ki2];
       normal[1] = ((double *) normals -> data)[ki2+1];
+
       //Edge flux computation
       flux_function(ql, qr, zl, zr,
@@ -988,37 +624,36 @@
 		    epsilon, g,
 		    edgeflux, &max_speed);
-      //update triangle k
-      ((long *) already_computed_flux->data)[ki]=call;
-      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
-      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
-      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
-      //update the neighbour n
-      if (n>=0){
-	((long *) already_computed_flux->data)[nm]=call;
-	((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
-	((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
-	((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
+
+
+      //flux -= edgeflux * edgelengths[k,i]
+      for (j=0; j<3; j++) {
+	flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
       }
-      ///for (j=0; j<3; j++) {
-	///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
-	///}
-	//Update timestep
-	//timestep = min(timestep, domain.radii[k]/max_speed)
-	//FIXME: SR Add parameter for CFL condition
-	if (max_speed > epsilon) {
-	  timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
-	  //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
-	  if (n>=0)
-	    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
-	}
+
+      //Update timestep
+      //timestep = min(timestep, domain.radii[k]/max_speed)
+      //FIXME: SR Add parameter for CFL condition
+      if (max_speed > epsilon) {
+	timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
+      }
     } // end for i
+
     //Normalise by area and store for when all conserved
     //quantities get updated
-    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
-    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
-    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
+    // flux /= areas[k]
+    for (j=0; j<3; j++) {
+      flux[j] /= ((double *) areas -> data)[k];
+    }
+
+    ((double *) stage_explicit_update -> data)[k] = flux[0];
+    ((double *) xmom_explicit_update -> data)[k] = flux[1];
+    ((double *) ymom_explicit_update -> data)[k] = flux[2];
+
   } //end for k
+
   return Py_BuildValue("d", timestep);
 }
+
+
 
 PyObject *protect(PyObject *self, PyObject *args) {
@@ -1178,68 +813,6 @@
 }
 
-PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
-  //a faster version of h_limiter above
-  PyObject *domain, *Tmp;
-  PyArrayObject
-    *hv, *hc, //Depth at vertices and centroids
-    *hvbar,   //Limited depth at vertices (return values)
-    *neighbours;
-
-  int k, i, N, k3,k0,k1,k2;
-  int dimensions[2];
-  double beta_h; //Safety factor (see config.py)
-  double hmin, hmax, dh[3];
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
-    return NULL;
-  neighbours = get_consecutive_array(domain, "neighbours");
-
-  //Get safety factor beta_h
-  Tmp = PyObject_GetAttrString(domain, "beta_h");
-  if (!Tmp)
-    return NULL;
-  beta_h = PyFloat_AsDouble(Tmp);
-
-  Py_DECREF(Tmp);
-
-  N = hc -> dimensions[0];
-
-  //Create hvbar
-  dimensions[0] = N;
-  dimensions[1] = 3;
-  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
-  for (k=0;k<N;k++){
-    k3=k*3;
-    //get the ids of the neighbours
-    k0 = ((long*) neighbours -> data)[k3];
-    k1 = ((long*) neighbours -> data)[k3+1];
-    k2 = ((long*) neighbours -> data)[k3+2];
-    //set hvbar provisionally
-    for (i=0;i<3;i++){
-      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
-      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
-    }
-    hmin=((double*) hc -> data)[k];
-    hmax=hmin;
-    if (k0>=0){
-      hmin=min(hmin,((double*) hc -> data)[k0]);
-      hmax=max(hmax,((double*) hc -> data)[k0]);
-    }
-    if (k1>=0){
-      hmin=min(hmin,((double*) hc -> data)[k1]);
-      hmax=max(hmax,((double*) hc -> data)[k1]);
-    }
-    if (k2>=0){
-      hmin=min(hmin,((double*) hc -> data)[k2]);
-      hmax=max(hmax,((double*) hc -> data)[k2]);
-    }
-    hmin-=((double*) hc -> data)[k];
-    hmax-=((double*) hc -> data)[k];
-    limit_gradient(dh,hmin,hmax,beta_h);
-    for (i=0;i<3;i++)
-      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
-  }
-  return PyArray_Return(hvbar);
-}
+
+
 
 PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
@@ -1292,5 +865,4 @@
 
   {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
-  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
   {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
   {"gravity", gravity, METH_VARARGS, "Print out"},
@@ -1299,6 +871,4 @@
    METH_VARARGS, "Print out"},
   {"h_limiter", h_limiter,
-   METH_VARARGS, "Print out"},
-  {"h_limiter_sw", h_limiter_sw,
    METH_VARARGS, "Print out"},
   {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
