Index: anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 4217)
+++ anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 4218)
@@ -147,4 +147,5 @@
 
         self.minimum_allowed_height = minimum_allowed_height
+        self.H0 = minimum_allowed_height  # Minimal height for flux computation        
         self.maximum_allowed_speed = maximum_allowed_speed
         self.g = g
@@ -832,5 +833,8 @@
          compute_fluxes_ext_central as compute_fluxes_ext
     
-    domain.timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.g,
+
+    domain.timestep = compute_fluxes_ext(timestep, domain.epsilon,
+                                         domain.H0,                                         
+                                         domain.g,
                                          domain.neighbours,
                                          domain.neighbour_edges,
@@ -2082,5 +2086,4 @@
                     protect_against_infinitesimal_and_negative_heights_c
 
-
     #distribute_to_vertices_and_edges =\
     #              distribute_to_vertices_and_edges_c #(like MH's)
Index: anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 4217)
+++ anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 4218)
@@ -129,5 +129,5 @@
 		  double z_left, double z_right,
 		  double n1, double n2,
-		  double epsilon, double g,
+		  double epsilon, double H0, double g,
 		  double *edgeflux, double *max_speed) {
 
@@ -152,4 +152,6 @@
   double flux_right[3], flux_left[3];
 
+  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
+  
   //Copy conserved quantities to protect from modification
   for (i=0; i<3; i++) {
@@ -173,5 +175,5 @@
     u_left = 0.0;
   } else {
-    u_left = uh_left/h_left;
+    u_left = uh_left/(h_left + h0/h_left);
   }
 
@@ -184,5 +186,5 @@
     u_right = 0.0;
   } else {
-    u_right = uh_right/h_right;
+    u_right = uh_right/(h_right + h0/h_right);
   }
 
@@ -249,8 +251,9 @@
 
 // Computational function for flux computation (using stage w=z+h)
+// FIXME (Ole): Is this used anywhere??
 int flux_function_kinetic(double *q_left, double *q_right,
 		  double z_left, double z_right,
 		  double n1, double n2,
-		  double epsilon, double g,
+		  double epsilon, double H0, double g,
 		  double *edgeflux, double *max_speed) {
 
@@ -270,4 +273,5 @@
   double q_left_copy[3], q_right_copy[3];
 
+  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
 
   //Copy conserved quantities to protect from modification
@@ -292,5 +296,5 @@
     u_left = 0.0;
   } else {
-    u_left = uh_left/h_left;
+    u_left = uh_left/(h_left + h0/h_left);
   }
 
@@ -303,5 +307,5 @@
     u_right = 0.0;
   } else {
-    u_right = uh_right/h_right;
+    u_right = uh_right/(h_right + h0/h_right);
   }
 
@@ -1175,4 +1179,5 @@
     domain.timestep = compute_fluxes(timestep,
                                      domain.epsilon,
+				     domain.H0,
                                      domain.g,
                                      domain.neighbours,
@@ -1221,5 +1226,5 @@
 
   //Local variables
-  double timestep, max_speed, epsilon, g;
+  double timestep, max_speed, epsilon, g, H0;
   double normal[2], ql[3], qr[3], zl, zr;
   double edgeflux[3]; //Work arrays for summing up fluxes
@@ -1231,7 +1236,8 @@
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOO",
 			&timestep,
 			&epsilon,
+			&H0,
 			&g,
 			&neighbours,
@@ -1300,5 +1306,5 @@
       flux_function_central(ql, qr, zl, zr,
 		    normal[0], normal[1],
-		    epsilon, g,
+		    epsilon, H0, g,
 		    edgeflux, &max_speed);
       //update triangle k
@@ -1363,4 +1369,5 @@
     domain.timestep = compute_fluxes(timestep,
                                      domain.epsilon,
+                                     domain.H0,
                                      domain.g,
                                      domain.neighbours,
@@ -1408,5 +1415,5 @@
 
   //Local variables
-  double timestep, max_speed, epsilon, g;
+  double timestep, max_speed, epsilon, g, H0;
   double normal[2], ql[3], qr[3], zl, zr;
   double edgeflux[3]; //Work arrays for summing up fluxes
@@ -1418,7 +1425,8 @@
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
 			&timestep,
 			&epsilon,
+			&H0,
 			&g,
 			&neighbours,
@@ -1485,5 +1493,5 @@
       flux_function_kinetic(ql, qr, zl, zr,
 		    normal[0], normal[1],
-		    epsilon, g,
+		    epsilon, H0, g,
 		    edgeflux, &max_speed);
       //update triangle k
Index: anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 4217)
+++ anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 4218)
@@ -103,4 +103,6 @@
             raise 'Should have raised an exception'
 
+
+    #FIXME (Ole): Individual flux tests do NOT test C implementation directly.    
     def test_flux_zero_case(self):
         ql = zeros( 3, Float )
@@ -594,4 +596,5 @@
 
 
+    # FIXME (Ole): Need test like this for fluxes in very shallow water.    
     def test_compute_fluxes3(self):
         #Random values, incl momentum
@@ -969,5 +972,6 @@
         domain = Domain(points, vertices, boundary) # Create domain
         domain.set_quantities_to_be_stored(None)
-        domain.set_maximum_allowed_speed(100) # 
+        domain.set_maximum_allowed_speed(100) #
+        domain.H0 = 0 # Backwards compatibility (6/2/7)
         
 
@@ -1050,4 +1054,5 @@
 
 
+        
         assert allclose(gauge_values[0], G0)
         assert allclose(gauge_values[1], G1)
@@ -2353,4 +2358,5 @@
         domain.beta_vh     = 0.9
         domain.beta_vh_dry = 0.9
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
         
         # Boundary conditions
@@ -2453,4 +2459,5 @@
         domain.smooth = False
         domain.default_order=1
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         # Boundary conditions
@@ -2503,4 +2510,5 @@
         domain.beta_vh_dry = 0.9        
         #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         # Boundary conditions
@@ -2566,4 +2574,5 @@
         domain.beta_vh_dry = 0.9        
         domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         # Boundary conditions
@@ -2615,5 +2624,6 @@
         domain.beta_uh_dry = 0.9
         domain.beta_vh     = 0.9
-        domain.beta_vh_dry = 0.9        
+        domain.beta_vh_dry = 0.9
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         # Boundary conditions
@@ -2799,4 +2809,5 @@
         domain.default_order=1
         domain.beta_h = 0.0 #Use first order in h-limiter
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         #Bed-slope and friction
@@ -2960,4 +2971,5 @@
         domain.beta_vh_dry = 0.9
         domain.beta_h = 0.0 #Use first order in h-limiter
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -3058,4 +3070,5 @@
         domain.beta_vh_dry = 0.9
         domain.beta_h = 0.0 #Use first order in h-limiter
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -3152,4 +3165,5 @@
         domain.beta_vh_dry = 0.9
         domain.beta_h = 0.0 #Use first order in h-limiter
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -3255,4 +3269,134 @@
 
 
+
+    def test_bedslope_problem_second_order_more_steps_feb_2007(self):
+        """test_bedslope_problem_second_order_more_steps_feb_2007
+
+        Test shallow water finite volumes, using parameters from feb 2007 rather
+        than backward compatibility ad infinitum
+        
+        """
+        from mesh_factory import rectangular
+        from Numeric import array
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order = 2
+        domain.beta_w      = 0.9
+        domain.beta_w_dry  = 0.9
+        domain.beta_uh     = 0.9
+        domain.beta_uh_dry = 0.9
+        domain.beta_vh     = 0.9
+        domain.beta_vh_dry = 0.9
+        domain.beta_h = 0.0 #Use first order in h-limiter
+        domain.H0 = 0.001 
+
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        def x_slope(x, y):
+            return -x/3
+
+        domain.set_quantity('elevation', x_slope)
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        #Initial condition
+        domain.set_quantity('stage', expression = 'elevation + 0.05')
+        domain.check_integrity()
+
+        assert allclose(domain.quantities['stage'].centroid_values,
+                        [0.01296296, 0.03148148, 0.01296296,
+                        0.03148148, 0.01296296, 0.03148148,
+                        0.01296296, 0.03148148, 0.01296296,
+                        0.03148148, 0.01296296, 0.03148148,
+                        -0.04259259, -0.02407407, -0.04259259,
+                        -0.02407407, -0.04259259, -0.02407407,
+                        -0.04259259, -0.02407407, -0.04259259,
+                        -0.02407407, -0.04259259, -0.02407407,
+                        -0.09814815, -0.07962963, -0.09814815,
+                        -0.07962963, -0.09814815, -0.07962963,
+                        -0.09814815, -0.07962963, -0.09814815,
+                        -0.07962963, -0.09814815, -0.07962963,
+                        -0.1537037 , -0.13518519, -0.1537037,
+                        -0.13518519, -0.1537037, -0.13518519,
+                        -0.1537037 , -0.13518519, -0.1537037,
+                        -0.13518519, -0.1537037, -0.13518519,
+                        -0.20925926, -0.19074074, -0.20925926,
+                        -0.19074074, -0.20925926, -0.19074074,
+                        -0.20925926, -0.19074074, -0.20925926,
+                        -0.19074074, -0.20925926, -0.19074074,
+                        -0.26481481, -0.2462963, -0.26481481,
+                        -0.2462963, -0.26481481, -0.2462963,
+                        -0.26481481, -0.2462963, -0.26481481,
+                        -0.2462963, -0.26481481, -0.2462963])
+
+
+        #print domain.quantities['stage'].extrapolate_second_order()
+        #domain.distribute_to_vertices_and_edges()
+        #print domain.quantities['stage'].vertex_values[:,0]
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
+            pass
+
+
+        assert allclose(domain.quantities['stage'].centroid_values,
+     [-0.02907614, -0.01462948, -0.02971293, -0.01435568, -0.02930822, -0.0141715,
+      -0.02900256, -0.01402774, -0.02897007, -0.01407298, -0.02958336, -0.01393198,
+      -0.07599907, -0.06253965, -0.07666738, -0.06312924, -0.07639753, -0.06265574,
+      -0.07572658, -0.06235626, -0.07569544, -0.0624551,  -0.07653298, -0.06289883,
+      -0.1236842,  -0.11090096, -0.12238737, -0.11116539, -0.12190705, -0.11072785,
+      -0.1208281,  -0.11001521, -0.12039674, -0.11011275, -0.1210331,  -0.11013222,
+      -0.16910086, -0.1583269,  -0.16731364, -0.15787261, -0.16655826, -0.15698876,
+      -0.16497539, -0.15560667, -0.16339389, -0.15509626, -0.16364638, -0.15425332,
+      -0.18773952, -0.19904678, -0.18906095, -0.1985912,  -0.18703781, -0.19698407,
+      -0.18337965, -0.19506463, -0.18190047, -0.19418413, -0.18587491, -0.19576587,
+      -0.1398922,  -0.14172812, -0.14134412, -0.14563187, -0.14097668, -0.14375562,
+      -0.13787839, -0.14035428, -0.13594691, -0.13938278, -0.13597595, -0.14218274])
+                      
+
+        assert allclose(domain.quantities['xmomentum'].centroid_values,
+     [ 0.00831991,  0.00327584,  0.0073476,   0.0034367,   0.00767331,  0.00356348,
+       0.00791201,  0.00364529,  0.00783242,  0.00349728,  0.0069742,   0.0031689,
+       0.02165378,  0.01421696,  0.02016895,  0.01318091,  0.02036575,  0.01369801,
+       0.02105492,  0.01400313,  0.02074176,  0.01356261,  0.01887467,  0.01232703,
+       0.03774493,  0.02854758,  0.03688063,  0.02759182,  0.03731697,  0.02811662,
+       0.03871488,  0.02912926,  0.03880077,  0.02803529,  0.0354565,   0.02599893,
+       0.06319926,  0.04729427,  0.05761625,  0.04591713,  0.05789906,  0.0468986,
+       0.05985286,  0.04870638,  0.06169141,  0.04811799,  0.05656696,  0.04415568,
+       0.08488718,  0.07187458,  0.07834844,  0.06842993,  0.07985979,  0.06981954,
+       0.08200942,  0.07216429,  0.08378261,  0.07273359,  0.08040488,  0.06646477,
+       0.01631518,  0.04691654,  0.02066439,  0.04441294,  0.02115705,  0.04560776,
+       0.02160783,  0.04664204,  0.02174952,  0.0479616,   0.02281756,  0.05667927])
+
+
+        assert allclose(domain.quantities['ymomentum'].centroid_values,
+     [  1.49908296e-04,  -3.32118806e-04,  -1.55139120e-04,  -2.97772609e-04,
+       -9.57477241e-05,  -3.11011790e-04,  -1.58896911e-04,  -3.76997605e-04,
+       -1.97659519e-04,  -3.34831296e-04,   6.54935308e-05,  -8.43493883e-06,
+        5.05063334e-04,  -1.43305846e-04,  -6.76061638e-05,  -5.01728590e-04,
+       -8.39003270e-05,  -4.64804117e-04,  -1.95135035e-04,  -5.88384357e-04,
+       -2.69800068e-04,  -5.35718006e-04,   2.59588334e-04,   3.00642498e-05,
+        5.15520349e-04,   1.05711270e-04,   9.26087960e-05,  -3.71855339e-04,
+        1.16386690e-04,  -3.82345749e-04,  -1.61741416e-04,  -6.31090292e-04,
+       -4.74530925e-04,  -6.95229436e-04,   6.08967544e-05,   2.20974020e-04,
+       -6.30000309e-04,   2.42320158e-04,  -5.89290588e-04,  -7.03283441e-05,
+       -4.18421888e-04,   6.62703090e-05,  -7.68647846e-04,  -3.40294284e-04,
+       -1.67585557e-03,  -7.40500723e-04,  -1.60020305e-03,   5.62070746e-05,
+       -1.48807206e-03,  -1.84455791e-03,  -2.27365819e-03,  -1.67381169e-03,
+       -1.95607481e-03,  -1.47497442e-03,  -1.73851968e-03,  -1.85314716e-03,
+       -2.01489344e-03,  -2.17608206e-03,  -1.66072261e-03,  -1.15856505e-03,
+       -1.18717624e-03,  -2.94595857e-03,  -3.59685615e-03,  -5.13811671e-03,
+       -6.17481232e-03,  -5.98871894e-03,  -6.00593324e-03,  -5.01282532e-03,
+       -4.51124363e-03,  -3.06579417e-03,   6.07680580e-04,  -4.80786234e-04])
+
+        os.remove(domain.get_name() + '.sww')
+
+
     def test_temp_play(self):
 
@@ -3274,4 +3418,5 @@
         domain.beta_vh_dry = 0.9
         domain.beta_h = 0.0 #Use first order in h-limiter
+        domain.H0 = 0 # Backwards compatibility (6/2/7)        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
