Index: anuga_work/development/anuga_1d/comp_flux_ext.c
===================================================================
--- anuga_work/development/anuga_1d/comp_flux_ext.c	(revision 5564)
+++ anuga_work/development/anuga_1d/comp_flux_ext.c	(revision 5587)
@@ -21,8 +21,36 @@
 	
 
+
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+double _compute_speed(double *uh, 
+		      double *h, 
+		      double epsilon, 
+		      double h0) {
+  
+  double u;
+  
+  if (*h < epsilon) {
+    *h = 0.0;  //Could have been negative
+     u = 0.0;
+  } else {
+    u = *uh/(*h + h0/ *h);    
+  }
+  
+
+  // Adjust momentum to be consistent with speed
+  *uh = u * *h;
+  
+  return u;
+}
+
+
+
 //Innermost flux function (using w=z+h)
 int _flux_function(double *q_left, double *q_right,
         double z_left, double z_right,
-		double normals, double g, double epsilon,
+		   double normals, double g, double epsilon, double h0,
 		double *edgeflux, double *max_speed) {
 		
@@ -33,5 +61,5 @@
 		double s_max, s_min, denom;
 		
-		
+		//printf("h0 = %f \n",h0);
 		ql[0] = q_left[0];
 		ql[1] = q_left[1];
@@ -54,20 +82,12 @@
 		h_left = w_left-z;
 		uh_left = ql[1];
-		if (h_left < epsilon) {
-			u_left = 0.0; h_left = 0.0;
-		}
-		else {
-			u_left = uh_left/h_left;
-		}
-	
+
+		u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
+
 		w_right = qr[0];
 		h_right = w_right-z;
 		uh_right = qr[1];
-		if (h_right < epsilon) {
-			u_right = 0.0; h_right = 0.0; 
-		}
-		else {
-			u_right = uh_right/h_right;
-		}
+
+		u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
   
 		soundspeed_left = sqrt(g*h_left);
@@ -113,21 +133,22 @@
 // Computational function for flux computation
 double _compute_fluxes_ext(double timestep,
-		double epsilon,
-		double g,
-		long* neighbours,
-		long* neighbour_vertices,
-		double* normals,
-		double* areas,
-		double* stage_edge_values,
-		double* xmom_edge_values,
-		double* bed_edge_values,
-		double* stage_boundary_values,
-		double* xmom_boundary_values,
-		double* stage_explicit_update,
-		double* xmom_explicit_update,
-		int number_of_elements,
-		double* max_speed_array) {
-		
-		double flux[2], ql[2], qr[2], edgeflux[2];
+			   double epsilon,
+			   double g,
+			   double h0,
+			   long* neighbours,
+			   long* neighbour_vertices,
+			   double* normals,
+			   double* areas,
+			   double* stage_edge_values,
+			   double* xmom_edge_values,
+			   double* bed_edge_values,
+			   double* stage_boundary_values,
+			   double* xmom_boundary_values,
+			   double* stage_explicit_update,
+			   double* xmom_explicit_update,
+			   int number_of_elements,
+			   double* max_speed_array) {
+		
+                double flux[2], ql[2], qr[2], edgeflux[2];
 		double zl, zr, max_speed, normal;
 		int k, i, ki, n, m, nm=0;
@@ -159,5 +180,5 @@
 				
 				normal = normals[ki];
-				_flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
+				_flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
 				flux[0] -= edgeflux[0];
 				flux[1] -= edgeflux[1];
@@ -209,12 +230,13 @@
 	*max_speed_array;
     
-  double timestep, epsilon, g;
+  double timestep, epsilon, g, h0;
   int number_of_elements;
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOiO",
 			&timestep,
 			&epsilon,
 			&g,
+			&h0,
 			&neighbours,
 			&neighbour_vertices,
@@ -240,4 +262,5 @@
 				 epsilon,
 				 g,
+				 h0,
 				 (long*) neighbours -> data,
 				 (long*) neighbour_vertices -> data,
@@ -281,5 +304,5 @@
 	*max_speed_array;
     
-  double timestep, epsilon, g;
+  double timestep, epsilon, g, h0;
   int number_of_elements;
 
@@ -299,4 +322,5 @@
     epsilon           = get_python_double(domain,"epsilon");
     g                 = get_python_double(domain,"g");
+    h0                = get_python_double(domain,"h0");
  
     
@@ -327,19 +351,20 @@
     // the explicit update arrays 
     timestep = _compute_fluxes_ext(timestep,
-				 epsilon,
-				 g,
-				 (long*) neighbours -> data,
-				 (long*) neighbour_vertices -> data,
-				 (double*) normals -> data,
-				 (double*) areas -> data,
-				 (double*) stage_vertex_values -> data,
-				 (double*) xmom_vertex_values -> data,
-				 (double*) bed_vertex_values -> data,
-				 (double*) stage_boundary_values -> data,
-				 (double*) xmom_boundary_values -> data,
-				 (double*) stage_explicit_update -> data,
-				 (double*) xmom_explicit_update -> data,
-				 number_of_elements,
-				 (double*) max_speed_array -> data);
+				   epsilon,
+				   g,
+				   h0,
+				   (long*) neighbours -> data,
+				   (long*) neighbour_vertices -> data,
+				   (double*) normals -> data,
+				   (double*) areas -> data,
+				   (double*) stage_vertex_values -> data,
+				   (double*) xmom_vertex_values -> data,
+				   (double*) bed_vertex_values -> data,
+				   (double*) stage_boundary_values -> data,
+				   (double*) xmom_boundary_values -> data,
+				   (double*) stage_explicit_update -> data,
+				   (double*) xmom_explicit_update -> data,
+				   number_of_elements,
+				   (double*) max_speed_array -> data);
 
 
