Index: /anuga_core/source/anuga/examples/run_sw_rectangle.py
===================================================================
--- /anuga_core/source/anuga/examples/run_sw_rectangle.py	(revision 3696)
+++ /anuga_core/source/anuga/examples/run_sw_rectangle.py	(revision 3696)
@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+#########################################################
+#
+#  Example showing improved wet dry limiting.
+#
+#
+#  Authors: Steve Roberts
+# October
+#
+#
+#
+#########################################################
+
+
+from Numeric import array
+from anuga.shallow_water import Domain
+
+
+# mesh partition routines
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
+
+
+M = 30
+points, vertices, boundary = rectangular(M, M, len1 = 1.0, len2 = 1.0)
+domain = Domain(points, vertices, boundary)
+print 'number of triangles = ', domain.number_of_elements
+
+
+#---------------------------
+#Boundaries
+#---------------------------
+from anuga.shallow_water import Reflective_boundary
+
+R = Reflective_boundary(domain)
+
+
+domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
+domain.check_integrity()
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0, h0=0.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.y0 = y0
+        self.y1 = y1
+        self.h  = h
+        self.h0 = h0
+
+    def __call__(self, x, y):
+        return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))
+
+domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
+
+import time
+t0 = time.time()
+
+
+# Turn on the visualisation
+
+rect = [0.0, 0.0, 1.0, 1.0]
+domain.initialise_visualiser()
+
+yieldstep = 0.002
+finaltime = 0.05
+
+
+#===============================================================================
+#Old Limiter
+#===============================================================================
+
+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
+
+#Check that the boundary value gets propagated to all elements
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+
+print 'That took %.2f seconds' %(time.time()-t0)
+print 'Note the small timesteps and the irregular flow'
+raw_input('press return to continue')
+
+#===============================================================================
+#New Limiter
+#===============================================================================
+
+t0 = time.time()
+domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
+domain.set_quantity('xmomentum', 0.0)
+domain.set_quantity('ymomentum', 0.0)
+domain.time = 0.0
+domain.default_order = 2
+domain.minimum_allowed_height = 0.001
+domain.beta_w      = 1.0
+domain.beta_w_dry  = 0.2
+domain.beta_uh     = 1.0
+domain.beta_uh_dry = 0.2
+domain.beta_vh     = 1.0
+domain.beta_vh_dry = 0.2
+
+#Check that the boundary value gets propagated to all elements
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+
+print 'That took %.2f seconds' %(time.time()-t0)
+print 'Note more uniform and large timesteps'
+raw_input('press return to continue')
+
+#===============================================================================
+#First Order
+#===============================================================================
+
+t0 = time.time()
+domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
+domain.set_quantity('xmomentum', 0.0)
+domain.set_quantity('ymomentum', 0.0)
+domain.time = 0.0
+
+domain.default_order = 1
+
+
+#Check that the boundary value gets propagated to all elements
+for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+    domain.write_time()
+
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
Index: /anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 3695)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 3696)
@@ -630,5 +630,5 @@
     *ymom_vertex_values,
 	*elevation_vertex_values;
-  PyObject *domain, *Tmp_w, *Tmp_w_dry, *Tmp_uh, *Tmp_uh_dry, *Tmp_vh, *Tmp_vh_dry;
+  PyObject *domain, *Tmp;
   //Local variables
   double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
@@ -639,4 +639,5 @@
   double hc, h0, h1, h2;
   double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
+  double minimum_allowed_height;
   //provisional jumps from centroids to v'tices and safety factor re limiting
   //by which these jumps are limited
@@ -661,39 +662,45 @@
 
   //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
-  Tmp_w = PyObject_GetAttrString(domain, "beta_w");
-  if (!Tmp_w)
-    return NULL;
-  beta_w = PyFloat_AsDouble(Tmp_w);
-  Py_DECREF(Tmp_w);
+  Tmp = PyObject_GetAttrString(domain, "beta_w");
+  if (!Tmp)
+    return NULL;
+  beta_w = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
   
-  Tmp_w_dry = PyObject_GetAttrString(domain, "beta_w_dry");
-  if (!Tmp_w_dry)
-    return NULL;
-  beta_w_dry = PyFloat_AsDouble(Tmp_w_dry);
-  Py_DECREF(Tmp_w_dry);
+  Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
+  if (!Tmp)
+    return NULL;
+  beta_w_dry = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
   
-  Tmp_uh = PyObject_GetAttrString(domain, "beta_uh");
-  if (!Tmp_uh)
-    return NULL;
-  beta_uh = PyFloat_AsDouble(Tmp_uh);
-  Py_DECREF(Tmp_uh);
+  Tmp = PyObject_GetAttrString(domain, "beta_uh");
+  if (!Tmp)
+    return NULL;
+  beta_uh = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
   
-  Tmp_uh_dry = PyObject_GetAttrString(domain, "beta_uh_dry");
-  if (!Tmp_uh_dry)
-    return NULL;
-  beta_uh_dry = PyFloat_AsDouble(Tmp_uh_dry);
-  Py_DECREF(Tmp_uh_dry); 
-
-  Tmp_vh = PyObject_GetAttrString(domain, "beta_vh");
-  if (!Tmp_vh)
-    return NULL;
-  beta_vh = PyFloat_AsDouble(Tmp_vh);
-  Py_DECREF(Tmp_vh);
+  Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
+  if (!Tmp)
+    return NULL;
+  beta_uh_dry = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp); 
+
+  Tmp = PyObject_GetAttrString(domain, "beta_vh");
+  if (!Tmp)
+    return NULL;
+  beta_vh = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
   
-  Tmp_vh_dry = PyObject_GetAttrString(domain, "beta_vh_dry");
-  if (!Tmp_vh_dry)
-    return NULL;
-  beta_vh_dry = PyFloat_AsDouble(Tmp_vh_dry);
-  Py_DECREF(Tmp_vh_dry);
+  Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
+  if (!Tmp)
+    return NULL;
+  beta_vh_dry = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
+  
+  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
+  if (!Tmp)
+    return NULL;
+  minimum_allowed_height = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);  
   
   number_of_elements = stage_centroid_values -> dimensions[0];
@@ -779,5 +786,5 @@
 	  hmin = qmin;
 	  beta_tmp = beta_w;
-	  if (hmin<0.001)
+	  if (hmin<minimum_allowed_height)
 		beta_tmp = beta_w_dry;
       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
@@ -804,5 +811,5 @@
       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
 	  beta_tmp = beta_uh;
-	  if (hmin<0.001)
+	  if (hmin<minimum_allowed_height)
 		beta_tmp = beta_uh_dry;
       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
@@ -829,5 +836,5 @@
       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
 	  beta_tmp = beta_vh;
-	  if (hmin<0.001)
+	  if (hmin<minimum_allowed_height)
 		beta_tmp = beta_vh_dry;
       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
Index: /anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py
===================================================================
--- /anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py	(revision 3695)
+++ /anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py	(revision 3696)
@@ -117,4 +117,11 @@
 domain.beta_vh_dry = 0.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
+
 yieldstep = 0.005
 finaltime = 1.0
