Index: trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py
===================================================================
--- trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py	(revision 9056)
@@ -73,5 +73,5 @@
 
         self.beta_w=1.0
-        self.beta_w_dry=0.4
+        self.beta_w_dry=0.0
         self.beta_uh=1.0
         self.beta_uh_dry=0.0
Index: trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c
===================================================================
--- trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c	(revision 9053)
+++ trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c	(revision 9056)
@@ -242,6 +242,8 @@
       if(i==0) edgeflux[i] += (s_max*s_min)*(max(q_right_rotated[i],ze) - max(q_left_rotated[i],ze));
       //if(i==0) edgeflux[i] += (s_max*s_min)*(hre - hle);
-      if(i==1) edgeflux[i] += (s_max*s_min)*(uh_right - uh_left);
-      if(i==2) edgeflux[i] += (s_max*s_min)*(vh_right - vh_left);
+      if((i==1)) edgeflux[i] += (s_max*s_min)*(uh_right - uh_left);
+      //if((i==1)) edgeflux[i] += (s_max*s_min)*(uh_right/(h_right+1.0e-08)*hre - uh_left/(h_left+1.0e-08)*hle);
+      if((i==2)) edgeflux[i] += (s_max*s_min)*(vh_right - vh_left);
+      //if((i==2)) edgeflux[i] += (s_max*s_min)*(vh_right/(h_right+1.0e-08)*hre - vh_left/(h_left+1.0e-08)*hle);
 
       edgeflux[i] *= inverse_denominator;
@@ -538,7 +540,7 @@
     for(k=0; k<number_of_elements; k++){
         hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
-        stage_explicit_update[k]=0.;
-        xmom_explicit_update[k]=0.;
-        ymom_explicit_update[k]=0.;
+        //stage_explicit_update[k]=0.;
+        //xmom_explicit_update[k]=0.;
+        //ymom_explicit_update[k]=0.;
 
         for(i=0;i<3;i++){
@@ -751,5 +753,5 @@
   double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2, dpth,momnorm;
   double dqv[3], qmin, qmax, hmin, hmax, bedmax,bedmin, stagemin;
-  double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight, tmp;
+  double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight, tmp, v0, v1, v2, vmax,vmin, vc;
   double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vel_norm, l1, l2, a_tmp, b_tmp, c_tmp,d_tmp;
   
@@ -845,21 +847,21 @@
   // condition) set its momentum to zero too. Aim is to prevent local 'pits'
   // with a non-zero water surface gradient building up momentum over time
-  for (k=0; k<number_of_elements;k++){
-
-      k3=k*3;
-      k0 = surrogate_neighbours[k3];
-      k1 = surrogate_neighbours[k3 + 1];
-      k2 = surrogate_neighbours[k3 + 2];
-
-      if((height_centroid_values[k0] < minimum_allowed_height | k0==k) &
-         (height_centroid_values[k1] < minimum_allowed_height | k1==k) &
-         (height_centroid_values[k2] < minimum_allowed_height | k2==k)){
-              xmom_centroid_store[k] = 0.;
-              xmom_centroid_values[k] = 0.;
-              ymom_centroid_store[k] = 0.;
-              ymom_centroid_values[k] = 0.;
-
-      }
-  }
+  //for (k=0; k<number_of_elements;k++){
+  //
+  //    k3=k*3;
+  //    k0 = surrogate_neighbours[k3];
+  //    k1 = surrogate_neighbours[k3 + 1];
+  //    k2 = surrogate_neighbours[k3 + 2];
+  //
+  //    if((height_centroid_values[k0] < minimum_allowed_height | k0==k) &
+  //       (height_centroid_values[k1] < minimum_allowed_height | k1==k) &
+  //       (height_centroid_values[k2] < minimum_allowed_height | k2==k)){
+  //            xmom_centroid_store[k] = 0.;
+  //            xmom_centroid_values[k] = 0.;
+  //            ymom_centroid_store[k] = 0.;
+  //            ymom_centroid_values[k] = 0.;
+
+  //    }
+  //}
 
   // Begin extrapolation routine
@@ -1052,4 +1054,15 @@
       hmin = min(min(h0, min(h1, h2)), hc);
       hmax = max(max(h0, max(h1, h2)), hc);
+
+      //vc = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k]+
+      //          ymom_centroid_values[k]*ymom_centroid_values[k]);
+      //v0 = sqrt(xmom_centroid_values[k0]*xmom_centroid_values[k0]+
+      //          ymom_centroid_values[k0]*ymom_centroid_values[k0]);
+      //v1 = sqrt(xmom_centroid_values[k1]*xmom_centroid_values[k1]+
+      //          ymom_centroid_values[k1]*ymom_centroid_values[k1]);
+      //v2 = sqrt(xmom_centroid_values[k2]*xmom_centroid_values[k2]+
+      //          ymom_centroid_values[k2]*ymom_centroid_values[k2]);
+      //vmax=max(max(v0*h0,max(v1*h1,v2*h2)),vc*hc);
+      //vmin=min(min(v0*h0,min(v1*h1,v2*h2)),vc*hc);
 
       // Look for strong changes in cell depth as an indicator of near-wet-dry
@@ -1075,4 +1088,6 @@
       // Set hfactor to 0 smoothly as hmin --> minumum_allowed_height 
       hfactor=min( 1.2*max(hmin-minimum_allowed_height,0.)/(max(hmin,0.)+1.*minimum_allowed_height), hfactor);
+
+      //hfactor=min(1.5*max(vmin,0.)/(max(vmax,0.)+1.0e-06), hfactor);
       //hfactor=min( 1.2*max(hc-minimum_allowed_height,0.)/(max(hc,0.)+1.*minimum_allowed_height), hfactor);
 
@@ -1252,5 +1267,5 @@
 
       // Limit the gradient
-      limit_gradient(dqv, qmin, qmax, beta_tmp);
+       limit_gradient(dqv, qmin, qmax, beta_tmp);
       
       for (i=0;i<3;i++)
Index: trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py	(revision 9056)
@@ -88,4 +88,5 @@
 domain.set_name('channel_floodplain1') # Output name
 #domain.set_store_vertices_uniquely(True)
+domain.set_flow_algorithm('DE1')
 #domain.use_edge_limiter=False
 #domain.extrapolate_velocity_second_order=False
Index: trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 9056)
@@ -1,4 +1,5 @@
 #import util
-from bal_and import plot_utils as util
+#from bal_and import plot_utils as util
+from anuga.utilities import plot_utils as util
 from matplotlib import pyplot as pyplot
 #import pylab
@@ -7,4 +8,5 @@
 #index=40
 index=570
+#index=225
 
 #p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
@@ -13,5 +15,5 @@
 
 #p2 = util.get_output('channel_floodplain1_bal_dev_lowvisc.sww', 0.01)
-p2 = util.get_output('channel_floodplain1.sww', 0.01)
+p2 = util.get_output('channel_floodplain1.sww', 0.001)
 p=util.get_centroids(p2, velocity_extrapolation=True)
 
Index: trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 9056)
@@ -28,8 +28,9 @@
 #domain.set_store_centroids(True)
 #domain.set_timestepping_method('euler')
-domain.set_flow_algorithm('DE1')
+#domain.set_flow_algorithm('DE1')
+domain.set_flow_algorithm('tsunami')
 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 
          'ymomentum': 2, 'elevation': 2})
-domain.minimum_storable_height=1.0e-03
+#domain.minimum_storable_height=1.0e-03
 domain.set_store_vertices_uniquely()
 #------------------
Index: trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py	(revision 9056)
@@ -5,5 +5,6 @@
 #---------------
 #import util
-from bal_and import plot_utils as util
+#from bal_and import plot_utils as util
+from anuga.utilities import plot_utils as util
 from matplotlib import pyplot as pyplot
 import numpy
Index: trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_2plot.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_2plot.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_2plot.py	(revision 9056)
@@ -9,6 +9,6 @@
 import scipy
 from matplotlib import pyplot as pyplot
-#from anuga.utilities import plot_utils as util
-from bal_and import plot_utils as util
+from anuga.utilities import plot_utils as util
+#from bal_and import plot_utils as util
 #--------------
 # Get variables
@@ -61,5 +61,5 @@
     pyplot.ylabel('depth (m)')
 
-pyplot.title('Flow depth at 400s -- there should be a central region with steady uniform flow ' )
+pyplot.title('Flow depth at end of simulation -- there should be a central region with steady uniform flow ' )
 pyplot.savefig('final_depth_v2.png')
 
@@ -83,9 +83,10 @@
 v=(p2.x==p2.x[tmp])
 
+ll=p2.vel.shape[0]-1
 #--------------------------------------------------
 # Check if y velocity is zero everywhere (it should be)
 #--------------------------------------------------
 pyplot.clf()
-pyplot.plot(p2.y[v],p2.yvel[100,v],'o', label='Numerical')
+pyplot.plot(p2.y[v],p2.yvel[ll,v],'o', label='Numerical')
 pyplot.plot((0,100),(0.0, 0.0),label='Analytical')
 pyplot.xlabel('Distance along the line x==50 (across the slope) (m)')
@@ -97,5 +98,5 @@
 
 pyplot.clf()
-pyplot.plot(p2.y[v],p2.xvel[100,v], 'o', label='Numerical')
+pyplot.plot(p2.y[v],p2.xvel[ll,v], 'o', label='Numerical')
 pyplot.plot((0,100),(uana,uana),label='Analytical')
 #pyplot.ylim([1.0,3.0])
Index: trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py	(revision 9056)
@@ -37,5 +37,5 @@
 #------------------------------------------------------------------------------
 def topography(x, y):
-	return -x/10.  + 1.*(numpy.sin(x/10.) +abs(y-50.)/10.) -0.*(x>80.) # linear bed slope
+	return -x/10.  #+ 1.*(numpy.sin(x/10.) +abs(y-50.)/10.) -0.*(x>80.) # linear bed slope
 
 def stagetopo(x,y):
Index: trunk/anuga_work/development/gareth/tests/wave_runup/plotme.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/wave_runup/plotme.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/wave_runup/plotme.py	(revision 9056)
@@ -66,7 +66,11 @@
 for i in range(len(model_shore_x)):
     # Compute index of shoreline
-    vtmp = p2.stage[i,:]>p2.elev+shoreline_depth_thresh1
+    #vtmp = p2.stage[i,:]>p2.elev+shoreline_depth_thresh1
+    #model_shore_x[i]=p2.x[vtmp].min()
     
-    model_shore_x[i]=p2.x[vtmp].min()
+    vtmp = p2.stage[i,:]<p2.elev+shoreline_depth_thresh1
+    model_shore_x[i]=p2.x[vtmp].max()
+    vtmp2=(p2.x>model_shore_x[i])
+    model_shore_x[i]=p2.x[vtmp2].min()
     #v_ind=(p2.stage[i,vtmp]-p2.elev[vtmp]).argmin()
     # Store x and xvel
@@ -76,12 +80,17 @@
     # Extract shoreline velocity. It is tricky, to avoid getting
     # zero velocity points -- might be a better way to do it?
-    vtmp2= (p2.stage[i,:]>p2.elev+shoreline_depth_thresh1)*\
-           (p2.stage[i,:]<p2.elev+shoreline_depth_thresh2)
-    mloc=abs(p2.xvel[i,vtmp2]).argmax()
+    #vtmp2= (p2.stage[i,:]>p2.elev+shoreline_depth_thresh1)*\
+    #       (p2.stage[i,:]<p2.elev+shoreline_depth_thresh2)
+    #mloc=abs(p2.xvel[i,vtmp2]).argmax()
+
     #print abs(p2.xvel[i,vtmp2]).max(), mloc
-    model_shore_u[i]=p2.xvel[i,vtmp2][mloc]
+    #model_shore_u[i]=p2.xvel[i,vtmp2][mloc]
     #vtmp2 = (p2.stage[i,:]>p2.elev+0.01)
     #mloc= (vtmp2*(p2.stage[i,:]-p2.elev)).argmin()
     #model_shore_u[i] = p2.xvel[i,mloc]
+    
+    vtmp2=vtmp
+    mloc=model_shore_x[i]
+    model_shore_u[i] = p2.xvel[i,p2.x==model_shore_x[i]].max()
 
 pyplot.plot(p2.time, model_shore_x-200.,'-o')
Index: trunk/anuga_work/development/gareth/tests/wave_runup/runup.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/wave_runup/runup.py	(revision 9053)
+++ trunk/anuga_work/development/gareth/tests/wave_runup/runup.py	(revision 9056)
@@ -10,6 +10,6 @@
 
 from math import sin, pi, exp
-from bal_and import *
-#from anuga import Domain as Domain
+#from bal_and import *
+from anuga import Domain as Domain
 
 #from anuga_tsunami import *
@@ -32,5 +32,5 @@
 domain.set_datadir('.')                          # Use current folder
 #domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
-#domain.set_flow_algorithm('DE1')
+domain.set_flow_algorithm('DE1')
 #------------------
 # Define topography
