Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1006)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1007)
@@ -696,11 +696,11 @@
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
+    hv = wv - zv
 
     #Call C-extension
     from shallow_water_ext import h_limiter
-    h_limiter(domain, hc, hv)
-
-    return hv
+    hvbar = h_limiter(domain, hc, hv)
+
+    return hvbar
 
 	    	                	    
Index: /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c	(revision 1006)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water_ext.c	(revision 1007)
@@ -239,4 +239,5 @@
       alpha = 1.0;  //Flat bed
 
+    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
 
     //	Let 
@@ -743,9 +744,9 @@
   PyArrayObject 
     *hv, *hc, //Depth at vertices and centroids        
-    //*hvbar,   //Limited depth at vertices (return values) FIXME: Not used?
+    *hvbar,   //Limited depth at vertices (return values) FIXME: Not used?
     *neighbours;
     
   int k, i, n, N, k3;
-  //int dimensions[2];
+  int dimensions[2];
   double beta_h; //Safety factor (see config.py)
   double *hmin, *hmax, hn;
@@ -769,7 +770,7 @@
   
   //Create hvbar
-  //dimensions[0] = N;
-  //dimensions[1] = 3;  
-  //hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+  dimensions[0] = N;
+  dimensions[1] = 3;  
+  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
  
 
@@ -785,5 +786,8 @@
     for (i=0; i<3; i++) {
       n = ((long*) neighbours -> data)[k3+i];
-      //((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k];  
+      
+      //Initialise hvbar with valus from hv
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];  
+      
       if (n >= 0) {
 	hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
@@ -796,5 +800,5 @@
   
   // Call underlying routine
-  _limit(N, beta_h, (double*) hc -> data, (double*) hv -> data, hmin, hmax);
+  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
 	 
   // // //Py_DECREF(domain); //FIXME: NEcessary?    	 
@@ -803,6 +807,6 @@
   
   //return result using PyArray to avoid memory leak  
-  //return PyArray_Return(hvbar);
-  return Py_BuildValue("");    
+  return PyArray_Return(hvbar);
+  //return Py_BuildValue("");    
 }
 
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1006)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 1007)
@@ -1445,5 +1445,81 @@
 	assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
 	assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
-	
+
+	
+	
+    def test_balance_deep_and_shallow(self):
+        """Test that balanced limiters preserve conserved quantites.
+	"""
+        import copy
+	
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        mesh = Domain(points, elements)        
+        mesh.check_integrity()
+
+        #Create a deliberate overshoot
+        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        mesh.set_quantity('elevation', 0) #Flat bed
+	stage = mesh.quantities['stage']     	
+	
+	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
+    
+        #Limit
+	mesh.distribute_to_vertices_and_edges()
+
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(mesh.number_of_elements):
+            assert allclose (ref_centroid_values[k],
+                             sum(stage.vertex_values[k,:])/3)
+        
+
+        #Now try with a non-flat bed - closely hugging initial stage in places
+	#This will create alphas in the range [0, 0.478260, 1]
+        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        mesh.set_quantity('elevation', [[0,0,0], 
+					[1.8,1.9,5.9], 
+					[4.6,0,0], 
+					[0,2,4]])
+	stage = mesh.quantities['stage']     	
+	
+	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
+	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy	
+    
+        #Limit
+	mesh.distribute_to_vertices_and_edges()
+
+	
+        #Assert that all vertex quantities have changed
+        for k in range(mesh.number_of_elements):
+	    #print ref_vertex_values[k,:], stage.vertex_values[k,:]	
+            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])		
+        #and assert that quantities are still conserved
+        from Numeric import sum
+        for k in range(mesh.number_of_elements):
+            assert allclose (ref_centroid_values[k],
+                             sum(stage.vertex_values[k,:])/3)
+        
+    
+        #Also check that Python and C version produce the same
+        assert allclose (stage.vertex_values,
+			 [[2,2,2],
+			  [1.93333333, 2.03333333, 6.03333333],
+			  [6.93333333, 4.53333333, 4.53333333],
+			  [5.33333333, 5.33333333, 5.33333333]])
+			  
+			 
+
+        
         
     def test_second_order_flat_bed_onestep(self):
@@ -2693,5 +2769,5 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(TestCase,'test')
-    #suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary_2')
+    #suite = unittest.makeSuite(TestCase,'test_balance_deep_and_shallow')
     runner = unittest.TextTestRunner()
     runner.run(suite)
