Index: trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py
===================================================================
--- trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py	(revision 8880)
+++ trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py	(revision 8884)
@@ -66,5 +66,5 @@
         self.set_CFL(1.0)
         self.set_use_kinematic_viscosity(False)
-        self.timestepping_method='rk2'#'rk2'#'rk2'#'euler'#'rk2' 
+        self.timestepping_method='rk2'#'euler'#'rk3'#'euler'#'rk2' 
         # The following allows storage of the negative depths associated with this method
         self.minimum_storable_height=-99999999999.0 
@@ -74,9 +74,9 @@
         self.extrapolate_velocity_second_order=True
 
-        self.beta_w=0.0
+        self.beta_w=1.0
         self.beta_w_dry=0.0
-        self.beta_uh=0.0
+        self.beta_uh=1.0
         self.beta_uh_dry=0.0
-        self.beta_vh=0.0
+        self.beta_vh=1.0
         self.beta_vh_dry=0.0
 
Index: trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c
===================================================================
--- trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c	(revision 8880)
+++ trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c	(revision 8884)
@@ -234,5 +234,6 @@
   // Flux formulas
   flux_left[0] = u_left*h_left;
-  //if(hc>h0){
+  //if(hc>h0 & hc_n > h0){
+  //if(w_left>z+h0 & w_right>z+h0){
     flux_left[1] = u_left*uh_left; //+ 0.5*g*h_left*h_left;
     flux_left[2] = u_left*vh_left;
@@ -243,5 +244,5 @@
 
   flux_right[0] = u_right*h_right;
-  //if(hc_n>h0){
+  //if(w_left>z+h0 & w_right>z+h0){
     flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right;
     flux_right[2] = u_right*vh_right;
@@ -253,5 +254,5 @@
   // Flux computation
   denom = s_max - s_min;
-  //if (denom < -epsilon) 
+  //if (denom < epsilon) 
   if (0==1) 
   { 
@@ -275,5 +276,5 @@
 
       // Standard smoothing term
-      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
+      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
 
       // Standard smoothing term, scaled
@@ -282,9 +283,9 @@
       // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
       // Physics 2:141-163
-      uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
-      t1 = (q_right_rotated[i] - uint);
-      t2 = (-q_left_rotated[i] + uint);
-      t3 = _minmod(t1,t2);
-      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3);
+      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
+      //t1 = (q_right_rotated[i] - uint);
+      //t2 = (-q_left_rotated[i] + uint);
+      //t3 = _minmod(t1,t2);
+      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3);
 
       // test this
@@ -457,19 +458,4 @@
             }
            
-            // GD HACK -- although I think this is a common technique 
-            //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
-            // unless the local centroid value is smaller 
-            //if(n>=0){
-            //    if((hc<=H0)&&(hc_n>H0)){
-            //        ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
-            //    }
-            //    if((hc_n<=H0)&&(hc>H0)){
-            //        qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
-            //    }
-
-            //}else{
-            //    // Treat the boundary case
-            //    // ??
-            //}
 
             smooth=1.0;
@@ -500,7 +486,13 @@
             // bedslope_work contains all gravity related terms -- weighting of
             // conservative and non-conservative versions.
+
+            if(ql[0]>zl){
+               tmp=1.0;
+            }else{
+               tmp=1.;
+            }
             if(hc> h0){
             //if(stage_centroid_values[k] > count_wet_edges[k]){
-              bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
+              bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl))*tmp + pressure_flux);
               //bedslope_work = length*(g*(hc*stage_centroid_values[k]*0.-
               //                           0.5*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)) 
@@ -509,10 +501,10 @@
               pressuregrad_store[ki] = bedslope_work;
             }else{
-              // NOTE: When hc<h0, pressure_flux is entirely computed from the
-              // other side of the edge. Hence, we can't necessarily satisfy
-              // well-balancing by just assuming that the 2nd and 3rd terms in
-              // bedslope_work cancel. So, we force it here -- equivalent to
-              // using the water surface gravity term
-              bedslope_work = length*(g*(hc*ql[0]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
+              // NOTE: When hc<h0, local contribution to the pressure_flux is
+              // entirely computed from the other side of the edge. Hence, we
+              // can't necessarily satisfy well-balancing by just assuming that
+              // the 2nd and 3rd terms in bedslope_work cancel. So, we force it
+              // here -- equivalent to using the water surface gravity term
+              bedslope_work = length*(g*(hc*ql[0]*tmp-0.0*max(ql[0]-zl,0.)*(ql[0]-zl))*tmp + 0.0*pressure_flux);
               //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
               //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
@@ -532,7 +524,12 @@
                 edgeflux_store[nm3 + 2 ] = edgeflux[2];
 
+                if(qr[0]>zr){
+                   tmp=1.0;
+                }else{
+                   tmp=1.;
+                }
                 if(hc_n> h0){
                 //if(stage_centroid_values[n] > count_wet_edges[n]){
-                  bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
+                  bedslope_work = length*(g*(hc_n*qr[0]*tmp-0.5*max(qr[0]-zr,0.)*(qr[0]-zr))*tmp + pressure_flux);
                   //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]*0.-0.5*max(stage_centroid_values[n]-zr,0.)*
                   //                                (stage_centroid_values[n]-zr)) + pressure_flux);
@@ -545,5 +542,5 @@
                   // bedslope_work cancel. So, we force it here -- equivalent to
                   // using the water surface gravity term
-                  bedslope_work = length*(g*(hc_n*qr[0]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
+                  bedslope_work = length*(g*(hc_n*qr[0]*tmp-0.0*max(qr[0]-zr,0.)*(qr[0]-zr))*tmp + 0.0*pressure_flux);
                   //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
                   //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
@@ -778,5 +775,5 @@
   // distance between the bed_centroid_value and the max bed_edge_value of
   // every triangle.
-  double minimum_relative_height=0.1; 
+  double minimum_relative_height=0.05; 
   int mass_added=0;
 
@@ -786,5 +783,6 @@
       hc = wc[k] - zc[k];
       // Definine the maximum bed edge value on triangle k.
-      bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+      //bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+      bmax = max(zv[3*k],max(zv[3*k+2],zv[3*k+1]));
 
       //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
@@ -796,7 +794,8 @@
       //if (hc < minimum_allowed_height*2.0 | hc < minimum_relative_height*(bmax-zc[k]) ){
       if (hc < minimum_allowed_height*2.0 ){
+      //if (hc < minimum_allowed_height*2.0 ){
             // Set momentum to zero and ensure h is non negative
-            xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
-            ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
+            //xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
+            //ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
             //xmomc[k] = xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
             //ymomc[k] = ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
@@ -809,5 +808,5 @@
              //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
              //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
-             bmin=zc[k];//-1.0e-03;
+             bmin=zc[k]-minimum_allowed_height;//-1.0e-03;
              // Minimum allowed stage = bmin
 
@@ -1006,4 +1005,21 @@
   }
 
+
+  // Alternative 'PROTECT' step
+  for(k=0; k<number_of_elements;k++){ 
+    if(count_wet_neighbours[k]<=3){
+        bedmax = max(elevation_centroid_values[k], 
+                     max(elevation_vertex_values[3*k],
+                         max(elevation_vertex_values[3*k+1], 
+                             elevation_vertex_values[3*k+2])));
+        if(cell_wetness_scale[k]==0. ||
+            stage_centroid_values[k] - elevation_centroid_values[k] < 0.05*(bedmax-elevation_centroid_values[k])){
+            xmom_centroid_store[k]=0.;
+            xmom_centroid_values[k]=0.;
+            ymom_centroid_store[k]=0.;
+            ymom_centroid_values[k]=0.;
+        }
+    } 
+  }
   
   // First treat all 'fully wet' cells, then treat 'partially dry' cells
@@ -1128,12 +1144,12 @@
                            max(elevation_centroid_values[k1], 
                                elevation_centroid_values[k2])));
-          bedmin = min(elevation_centroid_values[k], 
-                       min(elevation_centroid_values[k0],
-                           min(elevation_centroid_values[k1], 
-                               elevation_centroid_values[k2])));
-          //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]), 
-          //               min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
-          //                   min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
-          //                       max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
+          //bedmin = min(elevation_centroid_values[k], 
+          //             min(elevation_centroid_values[k0],
+          //                 min(elevation_centroid_values[k1], 
+          //                     elevation_centroid_values[k2])));
+          stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]), 
+                         min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
+                             min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
+                                 max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
 
           // Get the auxiliary triangle's vertex coordinates 
@@ -1163,5 +1179,5 @@
            
           //// Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
-          if ((area2 <= 0.))// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
+          if ((area2 <= 0.) )// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
           {
 
@@ -1217,5 +1233,5 @@
 
             //}   
-              //if( (k==k0) & (k==k1) & (k==k2)){
+              if( (k==k0) & (k==k1) & (k==k2)){
               //    // Isolated wet cell  -- use constant depth extrapolation for stage
               //    stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
@@ -1228,30 +1244,30 @@
               //    //stage_edge_values[k3+2] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2]);
               //    // Isolated wet cell -- constant stage extrapolation 
-              //    stage_edge_values[k3]   = stage_centroid_values[k];
-              //    stage_edge_values[k3+1] = stage_centroid_values[k];
-              //    stage_edge_values[k3+2] = stage_centroid_values[k];
-              //}else{
+                  stage_edge_values[k3]   = stage_centroid_values[k];
+                  stage_edge_values[k3+1] = stage_centroid_values[k];
+                  stage_edge_values[k3+2] = stage_centroid_values[k];
+              }else{
               //    // Cell with one wet neighbour. 
               //    // Find which neighbour is wet [if k!=k0, then k0 is wet]
-              //    if(k!=k0){
-              //      ktmp=k0;
-              //      ii=0;
-              //      xtmp = x0;
-              //      ytmp = y0;
-              //    }else if(k!=k1){
-              //      ktmp=k1;
-              //      ii=1;
-              //      xtmp = x1;
-              //      ytmp = y1;
-              //    }else if(k!=k2){
-              //      ktmp=k2;
-              //      ii=2;
-              //      xtmp = x2;
-              //      ytmp = y2;
-              //    }else{
-              //        printf("ERROR in extrapolation routine: Should have had one ki == k \n");
-              //        report_python_error(AT, "Negative triangle area");
-              //        return -1;
-              //    }
+                  if(k!=k0){
+                    ktmp=k0;
+                    ii=0;
+                    xtmp = x0;
+                    ytmp = y0;
+                  }else if(k!=k1){
+                    ktmp=k1;
+                    ii=1;
+                    xtmp = x1;
+                    ytmp = y1;
+                  }else if(k!=k2){
+                    ktmp=k2;
+                    ii=2;
+                    xtmp = x2;
+                    ytmp = y2;
+                  }else{
+                      printf("ERROR in extrapolation routine: Should have had one ki == k \n");
+                      report_python_error(AT, "Negative triangle area");
+                      return -1;
+                  }
 
               //    //if(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
@@ -1275,7 +1291,7 @@
               //           // stage_edge_values[k3+1] = stage_centroid_values[ktmp];
               //           // stage_edge_values[k3+2] = stage_centroid_values[ktmp];
-              //        //    //stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
-              //        //    //stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
-              //        //    //stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
+                     stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
+                     stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
+                     stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
               //        //}
 
@@ -1297,9 +1313,15 @@
               //    }
               //
-              //}
+              }
               // First order momentum / velocity extrapolation
-              stage_edge_values[k3]    = stage_centroid_values[k];
-              stage_edge_values[k3+1]  = stage_centroid_values[k];
-              stage_edge_values[k3+2]  = stage_centroid_values[k];
+             //if(stagemin>=bedmax){
+             //   stage_edge_values[k3]    = stage_centroid_values[k];
+             //   stage_edge_values[k3+1]  = stage_centroid_values[k];
+             //   stage_edge_values[k3+2]  = stage_centroid_values[k];
+             //}else{
+             //   stage_edge_values[k3]= stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3];
+             //   stage_edge_values[k3+1]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1];
+             //   stage_edge_values[k3+2]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2];
+             //}
               xmom_edge_values[k3]    = xmom_centroid_values[k];
               xmom_edge_values[k3+1]  = xmom_centroid_values[k];
@@ -1327,6 +1349,6 @@
             //hfactor=hmin/(hmin + 0.004);
           //}
-          hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height*2.)), 1.0);
           //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height)), 1.0);
+          hfactor= min(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height*0.1), 1.0);
           
           //-----------------------------------
@@ -1392,5 +1414,11 @@
           //if( (area2>0) ){ 
           //if(count_wet_neighbours[k]>0){
-          limit_gradient(dqv, qmin, qmax, beta_tmp);
+          //if(bedmax<=stagemin){
+             limit_gradient(dqv, qmin, qmax, beta_tmp);
+          //}else{
+          //  dqv[0]=-elevation_centroid_values[k] + elevation_edge_values[k3+0];
+          //  dqv[1]=-elevation_centroid_values[k] + elevation_edge_values[k3+1];
+          //  dqv[2]=-elevation_centroid_values[k] + elevation_edge_values[k3+2];
+          //}
           //printf("%lf \n ", beta_tmp);
           //}
@@ -1488,4 +1516,5 @@
             // (cell_wetness_scale[k]==1.0)){
           //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
+          
             limit_gradient(dqv, qmin, qmax, beta_tmp);
           //}else{
Index: trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 8880)
+++ trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 8884)
@@ -14,5 +14,5 @@
 #from swb2_domain import *
 #from balanced_basic import *
-#from balanced_dev import *
+from balanced_dev import *
 #from anuga_tsunami import *
 
@@ -22,8 +22,8 @@
 points, vertices, boundary = anuga.rectangular_cross(20,20, len1=1., len2=1.)
 
-domain=anuga.Domain(points,vertices,boundary)    # Create Domain
+domain=Domain(points,vertices,boundary)    # Create Domain
 domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
 #domain.set_timestepping_method('euler')
-domain.set_flow_algorithm('tsunami')
+#domain.set_flow_algorithm('tsunami')
 #------------------
 # Define topography
@@ -32,9 +32,9 @@
 #### Pathological
 #scale_me=100.0
-#boundary_ws=-0.2
+#boundary_ws=-0.2#1999
 #init_ws=-0.2
 #bumpiness=50. # Higher = shorter wavelength oscillations in topography
-#tstep=0.002
-#lasttime=1.1
+#tstep=0.2
+#lasttime=20.1
 
 ### Sensible
@@ -60,4 +60,5 @@
 domain.set_quantity('friction',0.00)             # Constant friction
 
+print '#-#', domain.minimum_allowed_height
 #def frict_change(x,y):
 #	return 0.2*(x>0.5)+0.1*(x<=0.5)
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 8880)
+++ trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py	(revision 8884)
@@ -43,6 +43,6 @@
     return stg#*(stg>topo) + topo*(stg<=topo)
 
-line1=[ [10.,0.], [10., 100.] ]
-Qin=20.
+line1=[ [9.,0.], [9., 100.] ]
+Qin=2.
 Inlet_operator(domain, line1,Qin)
 
