Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5305)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 5306)
@@ -306,9 +306,11 @@
 
         if self.default_order == 1:
-            self.set_timestepping_method('euler')
+            pass
+            #self.set_timestepping_method('euler')
             #self.set_all_limiters(beta_euler)
 
         if self.default_order == 2:
-            self.set_timestepping_method('rk2')
+            pass
+            #self.set_timestepping_method('rk2')
             #self.set_all_limiters(beta_rk2)
         
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5305)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 5306)
@@ -1409,8 +1409,14 @@
         extrapolate_from_gradient(self)
         
-    def extrapolate_second_order_and_limit(self):
+    def extrapolate_second_order_and_limit_by_edge(self):
         #Call correct module function
         #(either from this module or C-extension)
-        extrapolate_second_order_and_limit(self)
+        extrapolate_second_order_and_limit_by_edge(self)
+
+
+    def extrapolate_second_order_and_limit_by_vertex(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_and_limit_by_vertex(self)
 
     def bound_vertices_below_by_constant(self, bound):
@@ -1468,5 +1474,6 @@
          limit_gradient_by_neighbour,\
          extrapolate_from_gradient,\
-         extrapolate_second_order_and_limit,\
+         extrapolate_second_order_and_limit_by_edge,\
+         extrapolate_second_order_and_limit_by_vertex,\
          bound_vertices_below_by_constant,\
          bound_vertices_below_by_quantity,\
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 5305)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 5306)
@@ -266,5 +266,8 @@
 		     double* vertex_values,
 		     double* edge_values,
-		     long*   neighbours) {
+		     long*   neighbours,
+		     double* x_gradient,
+		     double* y_gradient) {
+
 
 	int i, k, k2, k3, k6;
@@ -305,6 +308,9 @@
 		    phi = min( min(r*beta, 1.0), phi);    
 		}
+
+		//Update gradient, vertex and edge values using phi limiter
+		x_gradient[k] = x_gradient[k]*phi;
+		y_gradient[k] = y_gradient[k]*phi;
     
-		//Update vertex and edge values using phi limiter
 		vertex_values[k3+0] = qc + phi*dqa[0];
 		vertex_values[k3+1] = qc + phi*dqa[1];
@@ -1126,5 +1132,5 @@
 
 
-PyObject *extrapolate_second_order_and_limit(PyObject *self, PyObject *args) {
+PyObject *extrapolate_second_order_and_limit_by_edge(PyObject *self, PyObject *args) {
     /* Compute edge values using second order approximation and limit values 
        so that edge values are limited by the two corresponding centroid values
@@ -1252,4 +1258,130 @@
 
 
+PyObject *extrapolate_second_order_and_limit_by_vertex(PyObject *self, PyObject *args) {
+    /* Compute edge values using second order approximation and limit values 
+       so that edge values are limited by the two corresponding centroid values
+       
+       Python Call:
+       extrapolate_second_order_and_limit(domain,quantity,beta)
+    */
+
+    PyObject *quantity, *domain;
+	
+    PyArrayObject
+	*domain_centroids,              //Coordinates at centroids
+	*domain_vertex_coordinates,     //Coordinates at vertices
+	*domain_number_of_boundaries,   //Number of boundaries for each triangle
+	*domain_surrogate_neighbours,   //True neighbours or - if one missing - self
+	*domain_neighbours,             //True neighbours, or if negative a link to boundary
+
+	*quantity_centroid_values,      //Values at centroids	
+	*quantity_vertex_values,        //Values at vertices
+	*quantity_edge_values,          //Values at edges
+	*quantity_phi,                  //limiter phi values
+	*quantity_x_gradient,           //x gradient
+	*quantity_y_gradient;           //y gradient
+    
+
+  // Local variables
+  int ntri;
+  double beta;
+  int err;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "O",&quantity)) {
+      PyErr_SetString(PyExc_RuntimeError, 
+		      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");	
+      return NULL;
+  }
+
+  domain = PyObject_GetAttrString(quantity, "domain");
+  if (!domain) {
+      PyErr_SetString(PyExc_RuntimeError, 
+		      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");	
+      return NULL;
+  }
+
+	
+  // Get pertinent variables
+  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
+  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
+  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
+  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
+  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
+
+  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
+  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
+  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
+  quantity_phi                  = get_consecutive_array(quantity, "phi");
+  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
+  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
+
+  beta = get_python_double(quantity,"beta");
+
+  ntri = quantity_centroid_values -> dimensions[0];
+
+  err = _compute_gradients(ntri,
+			   (double*) domain_centroids -> data,
+			   (double*) quantity_centroid_values -> data,
+			   (long*)   domain_number_of_boundaries -> data,
+			   (long*)   domain_surrogate_neighbours -> data,
+			   (double*) quantity_x_gradient -> data,
+			   (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ext.c: Internal function _compute_gradient failed");
+      return NULL;
+  }
+
+
+  err = _extrapolate_from_gradient(ntri, 
+				   (double*) domain_centroids -> data,
+				   (double*) quantity_centroid_values -> data,
+				   (double*) domain_vertex_coordinates -> data,
+				   (double*) quantity_vertex_values -> data,
+				   (double*) quantity_edge_values -> data,
+				   (double*) quantity_x_gradient -> data,
+				   (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
+      return NULL;
+  }
+
+
+  err = _limit_vertices_by_all_neighbours(ntri, beta,
+				       (double*) quantity_centroid_values -> data,
+				       (double*) quantity_vertex_values -> data,
+				       (double*) quantity_edge_values -> data,
+				       (long*)   domain_neighbours -> data,
+				       (double*) quantity_x_gradient -> data,
+				       (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+			  "quantity_ext.c: Internal function _limit_vertices_by_all_neighbours failed");
+      return NULL;
+  }
+
+
+  // Release
+  Py_DECREF(domain_centroids);
+  Py_DECREF(domain_surrogate_neighbours);
+  Py_DECREF(domain_number_of_boundaries);
+  Py_DECREF(domain_vertex_coordinates);
+
+  Py_DECREF(quantity_centroid_values);
+  Py_DECREF(quantity_vertex_values);
+  Py_DECREF(quantity_edge_values);
+  Py_DECREF(quantity_phi);
+  Py_DECREF(quantity_x_gradient);
+  Py_DECREF(quantity_y_gradient);
+
+  return Py_BuildValue("");
+}
+
+
 
 PyObject *compute_gradients(PyObject *self, PyObject *args) {
@@ -1445,5 +1577,7 @@
 	    *centroid_values, //Conserved quantities at centroids
 	    *edge_values,     //Conserved quantities at edges
-	    *neighbours;
+	    *neighbours,
+	    *x_gradient,
+	    *y_gradient;
 
 	double beta_w; //Safety factor
@@ -1481,5 +1615,8 @@
 	vertex_values    = get_consecutive_array(quantity, "vertex_values");
 	edge_values      = get_consecutive_array(quantity, "edge_values");
-	beta_w           = PyFloat_AsDouble(Tmp);
+	x_gradient       = get_consecutive_array(quantity, "x_gradient");
+	y_gradient       = get_consecutive_array(quantity, "y_gradient");
+	beta_w           = get_python_double(domain,"beta_w");
+
 
 
@@ -1490,5 +1627,9 @@
 						(double*) vertex_values -> data,
 						(double*) edge_values -> data,
-						(long*)   neighbours -> data);
+						(long*)   neighbours -> data,
+						(double*) x_gradient -> data,
+						(double*) y_gradient -> data);
+
+
 	
 	if (err != 0) {
@@ -1504,4 +1645,6 @@
 	Py_DECREF(vertex_values);
 	Py_DECREF(edge_values);
+	Py_DECREF(x_gradient);
+	Py_DECREF(y_gradient);
 	Py_DECREF(Tmp);
 
@@ -1557,5 +1700,5 @@
 	beta_w           = get_python_double(domain,"beta_w");
 
-	Py_DECREF(domain);
+
 
 	N = centroid_values -> dimensions[0];
@@ -1581,4 +1724,6 @@
 	Py_DECREF(vertex_values);
 	Py_DECREF(edge_values);
+	Py_DECREF(x_gradient);
+	Py_DECREF(y_gradient);
 
 
@@ -1625,5 +1770,4 @@
 
 
-	Py_DECREF(domain);
 
 	N = centroid_values -> dimensions[0];
@@ -1921,5 +2065,7 @@
 	{"extrapolate_from_gradient", extrapolate_from_gradient,
 		METH_VARARGS, "Print out"},
-	{"extrapolate_second_order_and_limit", extrapolate_second_order_and_limit,
+	{"extrapolate_second_order_and_limit_by_edge", extrapolate_second_order_and_limit_by_edge,
+		METH_VARARGS, "Print out"},
+	{"extrapolate_second_order_and_limit_by_vertex", extrapolate_second_order_and_limit_by_vertex,
 		METH_VARARGS, "Print out"},
 	{"interpolate_from_vertices_to_edges",
Index: /anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5305)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5306)
@@ -401,23 +401,8 @@
     def distribute_to_vertices_and_edges(self):
         # Call correct module function
-        # (either from this module or C-extension)
         if self.use_edge_limiter:
-            protect_against_infinitesimal_and_negative_heights(self)
-            for name in self.conserved_quantities:
-                Q = self.quantities[name]
-                if self._order_ == 1:
-                    Q.extrapolate_first_order()
-                elif self._order_ == 2:
-                    Q.extrapolate_second_order_and_limit()
-                else:
-                    raise 'Unknown order'
-            balance_deep_and_shallow(self)
-
-            #Compute edge values by interpolation
-            for name in self.conserved_quantities:
-                Q = self.quantities[name]
-                Q.interpolate_from_vertices_to_edges()
+            distribute_using_edge_limiter(self)            
         else:
-            distribute_to_vertices_and_edges(self)
+            distribute_using_vertex_limiter(self)
 
 
@@ -751,5 +736,5 @@
 
 
-def distribute_to_vertices_and_edges(domain):
+def distribute_using_vertex_limiter(domain):
     """Distribution from centroids to vertices specific to the
     shallow water wave
@@ -801,6 +786,5 @@
                 Q.extrapolate_first_order()
             elif domain._order_ == 2:
-                Q.extrapolate_second_order()
-                Q.limit()
+                Q.extrapolate_second_order_and_limit_by_vertex()
             else:
                 raise 'Unknown order'
@@ -808,4 +792,48 @@
 
     # Take bed elevation into account when water heights are small
+    balance_deep_and_shallow(domain)
+
+    # Compute edge values by interpolation
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        Q.interpolate_from_vertices_to_edges()
+
+
+
+def distribute_using_edge_limiter(domain):
+    """Distribution from centroids to edges specific to the
+    shallow water wave
+    equation.
+
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes by taking a weighted average between shallow
+    and deep cases.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    # Remove very thin layers of water
+    protect_against_infinitesimal_and_negative_heights(domain)
+
+
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        if domain._order_ == 1:
+            Q.extrapolate_first_order()
+        elif domain._order_ == 2:
+            Q.extrapolate_second_order_and_limit_by_edge()
+        else:
+            raise 'Unknown order'
+
     balance_deep_and_shallow(domain)
 
Index: /anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5305)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5306)
@@ -1823,6 +1823,9 @@
 }
 
-
-PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+PyObject *compute_fluxes_ext_central_new(PyObject *self, PyObject *args) {
   /*Compute all fluxes and the timestep suitable for all volumes
     in domain.
@@ -1840,86 +1843,81 @@
 
     Python call:
-    domain.timestep = compute_fluxes(timestep,
-                                     domain.epsilon,
-				     domain.H0,
-                                     domain.g,
-                                     domain.neighbours,
-                                     domain.neighbour_edges,
-                                     domain.normals,
-                                     domain.edgelengths,
-                                     domain.radii,
-                                     domain.areas,
-                                     tri_full_flag,
-                                     Stage.edge_values,
-                                     Xmom.edge_values,
-                                     Ymom.edge_values,
-                                     Bed.edge_values,
-                                     Stage.boundary_values,
-                                     Xmom.boundary_values,
-                                     Ymom.boundary_values,
-                                     Stage.explicit_update,
-                                     Xmom.explicit_update,
-                                     Ymom.explicit_update,
-                                     already_computed_flux,
-				     optimise_dry_cells)				     
+    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, bed)
 
 
     Post conditions:
       domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
+      returns timestep which is the largest step satisfying all volumes.
 
 
   */
 
-
-  PyArrayObject *neighbours, *neighbour_edges,
-    *normals, *edgelengths, *radii, *areas,
-    *tri_full_flag,
-    *stage_edge_values,
-    *xmom_edge_values,
-    *ymom_edge_values,
-    *bed_edge_values,
-    *stage_boundary_values,
-    *xmom_boundary_values,
-    *ymom_boundary_values,
-    *stage_explicit_update,
-    *xmom_explicit_update,
-    *ymom_explicit_update,
-    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
-    *max_speed_array; //Keeps track of max speeds for each triangle
-
-    
-  double timestep, epsilon, H0, g;
-  int optimise_dry_cells;
-    
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
-			&timestep,
-			&epsilon,
-			&H0,
-			&g,
-			&neighbours,
-			&neighbour_edges,
-			&normals,
-			&edgelengths, &radii, &areas,
-			&tri_full_flag,
-			&stage_edge_values,
-			&xmom_edge_values,
-			&ymom_edge_values,
-			&bed_edge_values,
-			&stage_boundary_values,
-			&xmom_boundary_values,
-			&ymom_boundary_values,
-			&stage_explicit_update,
-			&xmom_explicit_update,
-			&ymom_explicit_update,
-			&already_computed_flux,
-			&max_speed_array,
-			&optimise_dry_cells)) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
-    return NULL;
-  }
-
- 
+    PyObject 
+	*domain,
+	*stage, 
+	*xmom, 
+	*ymom, 
+	*bed;
+
+    PyArrayObject 
+	*neighbours, 
+	*neighbour_edges,
+	*normals, 
+	*edgelengths, 
+	*radii, 
+	*areas,
+	*tri_full_flag,
+	*stage_edge_values,
+	*xmom_edge_values,
+	*ymom_edge_values,
+	*bed_edge_values,
+	*stage_boundary_values,
+	*xmom_boundary_values,
+	*ymom_boundary_values,
+	*stage_explicit_update,
+	*xmom_explicit_update,
+	*ymom_explicit_update,
+	*already_computed_flux, //Tracks whether the flux across an edge has already been computed
+	*max_speed_array; //Keeps track of max speeds for each triangle
+
+    
+    double timestep, epsilon, H0, g;
+    int optimise_dry_cells;
+    
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
+	PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+	return NULL;
+    }
+
+    epsilon           = get_python_double(domain,"epsilon");
+    H0                = get_python_double(domain,"H0");
+    g                 = get_python_double(domain,"g");
+    optimise_dry_cells = get_python_integer(domain,"optimse_dry_cells");
+    
+    neighbours        = get_consecutive_array(domain, "neighbours");
+    neighbour_edges   = get_consecutive_array(domain, "neighbour_edges"); 
+    normals           = get_consecutive_array(domain, "normals");
+    edgelengths       = get_consecutive_array(domain, "edge_lengths");    
+    radii             = get_consecutive_array(domain, "radii");    
+    areas             = get_consecutive_array(domain, "areas");    
+    tri_full_flag     = get_consecutive_array(domain, "normals");
+    already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
+    max_speed_array   = get_consecutive_array(domain, "max_speed");
+    
+    stage_edge_values = get_consecutive_array(stage, "edge_values");    
+    xmom_edge_values  = get_consecutive_array(xmom, "edge_values");    
+    ymom_edge_values  = get_consecutive_array(ymom, "edge_values"); 
+    bed_edge_values   = get_consecutive_array(bed, "edge_values");    
+
+    stage_boundary_values = get_consecutive_array(stage, "boundary_values");    
+    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");    
+    ymom_boundary_values  = get_consecutive_array(ymom, "boundary_values"); 
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");    
+    ymom_explicit_update  = get_consecutive_array(ymom, "explicit_update"); 
+
+
   int number_of_elements = stage_edge_values -> dimensions[0];
 
@@ -1951,8 +1949,168 @@
 				     (double*) max_speed_array -> data,
 				     optimise_dry_cells);
+
+  Py_DECREF(neighbours);
+  Py_DECREF(neighbour_edges);
+  Py_DECREF(normals);
+  Py_DECREF(edgelengths);
+  Py_DECREF(radii);
+  Py_DECREF(areas);
+  Py_DECREF(tri_full_flag);
+  Py_DECREF(already_computed_flux);
+  Py_DECREF(max_speed_array);
+  Py_DECREF(stage_edge_values);
+  Py_DECREF(xmom_edge_values);
+  Py_DECREF(ymom_edge_values);
+  Py_DECREF(bed_edge_values);
+  Py_DECREF(stage_boundary_values);
+  Py_DECREF(xmom_boundary_values);
+  Py_DECREF(ymom_boundary_values);
+  Py_DECREF(stage_explicit_update);
+  Py_DECREF(xmom_explicit_update);
+  Py_DECREF(ymom_explicit_update);
+
   
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
 }
+
+
+
+
+
+
+PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function_central"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Python call:
+    domain.timestep = compute_fluxes(timestep,
+                                     domain.epsilon,
+				     domain.H0,
+                                     domain.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,
+                                     domain.radii,
+                                     domain.areas,
+                                     tri_full_flag,
+                                     Stage.edge_values,
+                                     Xmom.edge_values,
+                                     Ymom.edge_values,
+                                     Bed.edge_values,
+                                     Stage.boundary_values,
+                                     Xmom.boundary_values,
+                                     Ymom.boundary_values,
+                                     Stage.explicit_update,
+                                     Xmom.explicit_update,
+                                     Ymom.explicit_update,
+                                     already_computed_flux,
+				     optimise_dry_cells)				     
+
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+
+
+  */
+
+
+  PyArrayObject *neighbours, *neighbour_edges,
+    *normals, *edgelengths, *radii, *areas,
+    *tri_full_flag,
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *bed_edge_values,
+    *stage_boundary_values,
+    *xmom_boundary_values,
+    *ymom_boundary_values,
+    *stage_explicit_update,
+    *xmom_explicit_update,
+    *ymom_explicit_update,
+    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
+    *max_speed_array; //Keeps track of max speeds for each triangle
+
+    
+  double timestep, epsilon, H0, g;
+  int optimise_dry_cells;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
+			&timestep,
+			&epsilon,
+			&H0,
+			&g,
+			&neighbours,
+			&neighbour_edges,
+			&normals,
+			&edgelengths, &radii, &areas,
+			&tri_full_flag,
+			&stage_edge_values,
+			&xmom_edge_values,
+			&ymom_edge_values,
+			&bed_edge_values,
+			&stage_boundary_values,
+			&xmom_boundary_values,
+			&ymom_boundary_values,
+			&stage_explicit_update,
+			&xmom_explicit_update,
+			&ymom_explicit_update,
+			&already_computed_flux,
+			&max_speed_array,
+			&optimise_dry_cells)) {
+    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    return NULL;
+  }
+
+ 
+  int number_of_elements = stage_edge_values -> dimensions[0];
+
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  timestep = _compute_fluxes_central(number_of_elements,
+				     timestep,
+				     epsilon,
+				     H0,
+				     g,
+				     (long*) neighbours -> data,
+				     (long*) neighbour_edges -> data,
+				     (double*) normals -> data,
+				     (double*) edgelengths -> data, 
+				     (double*) radii -> data, 
+				     (double*) areas -> data,
+				     (long*) tri_full_flag -> data,
+				     (double*) stage_edge_values -> data,
+				     (double*) xmom_edge_values -> data,
+				     (double*) ymom_edge_values -> data,
+				     (double*) bed_edge_values -> data,
+				     (double*) stage_boundary_values -> data,
+				     (double*) xmom_boundary_values -> data,
+				     (double*) ymom_boundary_values -> data,
+				     (double*) stage_explicit_update -> data,
+				     (double*) xmom_explicit_update -> data,
+				     (double*) ymom_explicit_update -> data,
+				     (long*) already_computed_flux -> data,
+				     (double*) max_speed_array -> data,
+				     optimise_dry_cells);
+  
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
 
 
@@ -2490,4 +2648,5 @@
   {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
   {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
+  {"compute_fluxes_ext_central_new", compute_fluxes_ext_central_new, METH_VARARGS, "Print out"},
   {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
   {"gravity", gravity, METH_VARARGS, "Print out"},
Index: /anuga_core/source/anuga/utilities/util_ext.h
===================================================================
--- /anuga_core/source/anuga/utilities/util_ext.h	(revision 5305)
+++ /anuga_core/source/anuga/utilities/util_ext.h	(revision 5306)
@@ -298,4 +298,23 @@
 }
 
+int get_python_integer(PyObject *O, char *name) {
+  PyObject *TObject;
+  int tmp;
+  
+
+  //Get double from attribute
+  TObject = PyObject_GetAttrString(O, name);
+  if (!TObject) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object");
+    return 0;
+  }  
+  
+  tmp = PyFloat_AsDouble(TObject);
+  
+  Py_DECREF(TObject);
+  
+  return tmp;
+}
+
 
 PyObject *get_python_object(PyObject *O, char *name) {
Index: /anuga_core/source/anuga_parallel/run_parallel_advection.py
===================================================================
--- /anuga_core/source/anuga_parallel/run_parallel_advection.py	(revision 5305)
+++ /anuga_core/source/anuga_parallel/run_parallel_advection.py	(revision 5306)
@@ -102,19 +102,7 @@
     t0 = time.time()
 
-# Start the evolve computions
-import hotshot
-profiler = hotshot.Profile("hotshot." + str(numprocs) + "." + str(myid) + ".prof")
-s = '''for t in domain.evolve(yieldstep = 0.1, finaltime = 2.0):
-  if myid == 0:
-    domain.write_time()
-'''
-    
-
-result = profiler.runctx(s, globals(), locals())
-profiler.close()
-
-# for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
-#     if myid == 0:
-#         domain.write_time()
+for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+    if myid == 0:
+        domain.write_time()
 
 # Output some computation statistics
Index: /anuga_validation/convergence_study/plot_timeseries.py
===================================================================
--- /anuga_validation/convergence_study/plot_timeseries.py	(revision 5306)
+++ /anuga_validation/convergence_study/plot_timeseries.py	(revision 5306)
@@ -0,0 +1,141 @@
+## def animatesww2d_alt(sww_filename, movie_filename, range):
+##     """Plot cross section of model output
+##     """
+
+##     # Read SWW file	
+##     from Scientific.IO.NetCDF import NetCDFFile
+##     fid = NetCDFFile(sww_filename, 'r')
+    
+##     x = fid.variables['x']
+##     y = fid.variables['y']
+##     volumes = fid.variables['volumes']
+    
+##     elevation = fid.variables['elevation']
+##     time = fid.variables['time']
+##     stage = fid.variables['stage']
+##     xmomentum = fid.variables['xmomentum']
+##     ymomentum = fid.variables['ymomentum']
+    
+def animatesww2d(swwfile):
+    """Read in sww file and plot cross section of model output
+    """
+
+    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+
+    try:
+        fid = open(swwfile)
+    except Exception, e:
+        msg = 'File "%s" could not be opened: Error="%s"'\
+              %(swwfile, e)
+        raise msg
+
+    print 'swwfile', swwfile
+
+    # interpolation points are for y = 0
+    # number of increments in x
+    m = 1000 
+    max_x = 100000.
+    x = 0
+    points = []
+    x_points = []
+    for i in range(m):
+        x += max_x/m
+        points.append([x,0.])
+        x_points.append(x)
+            
+    time_thinning = 1
+    use_cache = True
+    verbose = True
+
+    from anuga.abstract_2d_finite_volumes.util import file_function
+    
+    f = file_function(swwfile,
+                      quantities = sww_quantity,
+                      interpolation_points = points,
+                      time_thinning = time_thinning,
+                      verbose = verbose,
+                      use_cache = use_cache)
+
+    n = len(f.get_time()) # number of time steps
+
+    from Numeric import zeros, Float
+    from math import sqrt
+    stage = zeros((n,m), Float)
+    elevation = zeros((n,m), Float)
+    xmom = zeros((n,m), Float)
+    ymom = zeros((n,m), Float)
+    momenta = zeros((n,m), Float)
+    speed = zeros((n,m), Float)
+    max_stages = zeros((n,m), Float)
+    max_stage = None
+    min_stages = zeros((n,m), Float)
+    min_stage = 100.
+    fid = open('stage.csv', 'w')
+    for k, location in enumerate(x_points):
+        for i, t in enumerate(f.get_time()):
+            w = f(t, point_id = k)[0]
+            z = f(t, point_id = k)[1]
+            uh = f(t, point_id = k)[2]
+            vh = f(t, point_id = k)[3]
+            depth = w-z      
+            m = sqrt(uh*uh + vh*vh)
+            if depth < 0.001:
+                vel = 0.0
+            else:
+                vel = m / (depth + 1.e-6/depth)
+
+            stage[i,k] = w
+            s = '%.2f,' %(w)
+            fid.write(s)
+            elevation[i,k] = z
+            xmom[i,k] = uh
+            ymom[i,k] = vh
+            momenta[i,k] = m
+            speed[i,k] = vel
+
+            if w > max_stage: max_stage = w
+            if w < min_stage: min_stage = w
+            max_stages[i,k] = max_stage
+            min_stages[i,k] = min_stage
+        s = '\n'
+        fid.write(s)
+
+    min_stage = min(min(min_stages))
+    max_stage = max(max(max_stages))
+    stage_axis = [0, max_x, min_stage-1, max_stage+1]
+    from pylab import plot, xlabel, ylabel, savefig, close, hold, axis, title,show
+    hold(False)
+    figs = []
+    j = 1
+    #max_vec = zeros(100, Float)
+
+    #if sys.platform == 'win32':
+    try:
+        for i in range(0,n,20):
+            #max_vec[i] = max(stage[i,:])
+            plot(x_points,stage[i,:]) #x_points,max_vec,'+') 
+            xlabel('x (m)')
+            ylabel('stage (m)')
+            axis(stage_axis)
+            name = 'time_%i' %i
+            savefig(name)
+            title(name)
+            figs.append(name)
+            show()
+    except:
+        pass
+    
+
+    close('all')
+
+if __name__ == '__main__':
+    import sys
+    swwfile = sys.argv[1]
+    
+    animatesww2d(swwfile)
+    
+#swwfile = join(getenv('INUNDATIONHOME'),'data','validation','convergence_study','5000','myexample2.sww')
+#animatesww2d(swwfile)
+
+# once png files have been created, they can then be dragged into
+# Windows Movie Maker and exported to a mpeg.
Index: /anuga_validation/convergence_study/wave.py
===================================================================
--- /anuga_validation/convergence_study/wave.py	(revision 5305)
+++ /anuga_validation/convergence_study/wave.py	(revision 5306)
@@ -39,4 +39,6 @@
 #start_screen_catcher(output_dir+sep)
 
+interactive_visualisation = False
+
 #------------------------------------------------------------------------------
 # Setup domain
@@ -77,5 +79,4 @@
 domain.beta_h    = 0.0
 
-interactive_visualisation = True
 
 
