Index: /anuga_core/documentation/user_manual/demos/channel3.py
===================================================================
--- /anuga_core/documentation/user_manual/demos/channel3.py	(revision 5314)
+++ /anuga_core/documentation/user_manual/demos/channel3.py	(revision 5315)
@@ -74,5 +74,5 @@
 # Evolve system through time
 #------------------------------------------------------------------------------
-for t in domain.evolve(yieldstep = 0.2, finaltime = 16.0):
+for t in domain.evolve(yieldstep = 0.1, finaltime = 16.0):
     print domain.timestepping_statistics()
 
Index: /anuga_core/source/anuga/config.py
===================================================================
--- /anuga_core/source/anuga/config.py	(revision 5314)
+++ /anuga_core/source/anuga/config.py	(revision 5315)
@@ -88,7 +88,9 @@
 tight_slope_limiters = True
 
-# Use centroid velocities to reconstruct momentum at vertices in very shallow water
+# Use centroid velocities to reconstruct momentum at vertices in
+# very shallow water
 # This option has a first order flavour to it, but we still have second order
-# reconstruction of stage and this option only applies in balance_deep_and_shallow when
+# reconstruction of stage and this option only applies in
+# balance_deep_and_shallow when
 # alpha < 1 so in deeper water the full second order scheme is used.
 #
Index: /anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5314)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 5315)
@@ -732,6 +732,5 @@
               Ymom.vertex_values,
               Elevation.vertex_values,
-              int(domain.optimise_dry_cells),
-              int(domain.use_centroid_velocities))
+              int(domain.optimise_dry_cells))
 
 
Index: /anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5314)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 5315)
@@ -550,13 +550,12 @@
 			      double H0,
 			      int tight_slope_limiters,
-			      int use_centroid_velocities,		      
+			      int use_centroid_velocities,
 			      double alpha_balance) {
 
-  int k, k3, i; //, excessive_froude_number=0;
+  int k, k3, i; 
 
   double dz, hmin, alpha, h_diff, hc_k;
   double epsilon = 1.0e-6; // FIXME: Temporary measure
   double hv[3]; // Depths at vertices
-  //double Fx, Fy; // Froude numbers
   double uc, vc; // Centroid speeds
 
@@ -651,7 +650,6 @@
 	alpha = 1.0;       
       }
-
     }
-	    
+
 
     //	Let
@@ -684,5 +682,25 @@
 
 	// Update momentum at vertices
-	if (use_centroid_velocities == 0) {
+	if (use_centroid_velocities == 1) {
+	  // This is a simple, efficient and robust option
+	  // It uses first order approximation of velocities, but retains
+	  // the order used by stage.
+	
+	  // Speeds at centroids
+	  if (hc_k > epsilon) {
+	    uc = xmomc[k]/hc_k;
+	    vc = ymomc[k]/hc_k;
+	  } else {
+	    uc = 0.0;
+	    vc = 0.0;
+	  }
+	  
+	  // Vertex momenta guaranteed to be consistent with depth guaranteeing
+	  // controlled speed
+	  hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
+	  xmomv[k3+i] = uc*hv[i];
+	  ymomv[k3+i] = vc*hv[i];
+	  
+	} else {
 	  // Update momentum as a linear combination of
 	  // xmomc and ymomc (shallow) and momentum
@@ -691,4 +709,8 @@
 	  // been established e.g. by the gradient limiter.
 
+	  // FIXME (Ole): I think this should be used with vertex momenta
+	  // computed above using centroid_velocities instead of xmomc 
+	  // and ymomc as they'll be more representative first order
+	  // values.
 	  
 	  xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
@@ -697,34 +719,9 @@
 	}
       }
-    } // If alpha == 1 use quantities as calculated by the gradient-limiter
-    
-    if (use_centroid_velocities == 1) {    
-      // This is a simple, efficient and robust option in shallow water
-      // It uses first order approximation of velocities, but retains
-      // the order used by stage.
-      // In this case the xmomv and ymomv calculations should be switched off
-      // in the gradient limiter.
-
-      for (i=0; i<3; i++) {    
-	
-	// Speeds at centroids
-	if (hc_k > epsilon) {
-	  uc = xmomc[k]/hc_k;
-	  vc = ymomc[k]/hc_k;
-	} else {
-	  uc = 0.0;
-	  vc = 0.0;
-	}
-	
-	// Vertex momenta guaranteed to be consistent with depth guaranteeing
-	// controlled speed
-	hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
-	xmomv[k3+i] = uc*hv[i];
-	ymomv[k3+i] = vc*hv[i];
-      } 
     }
   }
   return 0;
 }
+  
 
 
@@ -1020,6 +1017,5 @@
 				  double* ymom_vertex_values,
 				  double* elevation_vertex_values,
-				  int optimise_dry_cells,
-				  int use_centroid_velocities) {
+				  int optimise_dry_cells) {
 				  
 				  
@@ -1027,5 +1023,5 @@
   // Local variables
   double a, b; // Gradient vector used to calculate vertex values from centroids
-  int k,k0,k1,k2,k3,k6,coord_index, i;
+  int 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;
@@ -1192,16 +1188,4 @@
       for (i=0;i<3;i++)
 	stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
-      
-      
-      if (use_centroid_velocities == 1) {
-	// Use first order reconstruction using speeds only.
-	
-	// This happens in balance_deep_and_shallow so there
-	// is no need to do more here
-	
-	// FIXME (Ole): Optionally we could put the computation here but then
-	// it'd go into all other variants of the gradient limiter  
-	continue;
-      }
       
       
@@ -1500,10 +1484,10 @@
   double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;    
   double minimum_allowed_height, epsilon;
-  int optimise_dry_cells, number_of_elements, e, use_centroid_velocities;
+  int optimise_dry_cells, number_of_elements, e;
   
   // 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, "OOOOOOOOOOOOOii",
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
 			&domain,
 			&surrogate_neighbours,
@@ -1519,6 +1503,5 @@
 			&ymom_vertex_values,
 			&elevation_vertex_values,
-			&optimise_dry_cells,
-			&use_centroid_velocities)) {			
+			&optimise_dry_cells)) {			
 			
     PyErr_SetString(PyExc_RuntimeError, 
@@ -1565,6 +1548,5 @@
 				   (double*) ymom_vertex_values -> data,
 				   (double*) elevation_vertex_values -> data,
-				   optimise_dry_cells,
-				   use_centroid_velocities);
+				   optimise_dry_cells);
   if (e == -1) {
     // Use error string set inside computational routine
