Index: /trunk/anuga_core/source/anuga/structures/inlet_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/structures/inlet_operator.py	(revision 8548)
+++ /trunk/anuga_core/source/anuga/structures/inlet_operator.py	(revision 8549)
@@ -62,4 +62,6 @@
         volume = Q*timestep
         
+        #print Q, volume
+        
         assert current_volume + volume >= 0.0, 'Requesting too much water to be removed from an inlet!'
 
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 8548)
+++ /trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c	(revision 8549)
@@ -196,9 +196,9 @@
   // Flux formulas
   flux_left[0] = u_left*h_left;
-  flux_left[1] = u_left*uh_left + 0.0*g*h_left*h_left;
+  flux_left[1] = u_left*uh_left ; //+ 0.5*g*h_left*h_left;
   flux_left[2] = u_left*vh_left;
 
   flux_right[0] = u_right*h_right;
-  flux_right[1] = u_right*uh_right + 0.0*g*h_right*h_right;
+  flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right;
   flux_right[2] = u_right*vh_right;
 
@@ -221,5 +221,5 @@
       edgeflux[i] *= inverse_denominator;
     }
-
+    // Separate pressure flux, so we can apply different wet-dry hacks to it
     *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
     
@@ -288,5 +288,5 @@
 
     max_bed_edgevalue = malloc(number_of_elements*sizeof(double));
-    min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
+    //min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
     // Start computation
     call++; // Flag 'id' of flux calculation for this timestep
@@ -304,6 +304,6 @@
                                    max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
         //max_bed_edgevalue[k]= 0.5*(max_bed_edgevalue[j]+bed_centroid_values[k]);
-        min_bed_edgevalue[k] = min(bed_edge_values[3*k], 
-                                   min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
+        //min_bed_edgevalue[k] = min(bed_edge_values[3*k], 
+        //                           min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
     
     }
@@ -408,5 +408,5 @@
             // Prevent outflow from 'seriously' dry cells
             // Idea: The cell will not go dry if:
-            // mass_flux = edgeflux[0]*(dt*edgelength) <= Area_triangle*hc
+            // mass_flux = edgeflux[0]*(dt*edgelength) <= (1/3)Area_triangle*hc
             if((stage_centroid_values[k]<=max_bed_edgevalue[k])|
                (ql[0]<=zl)){
@@ -469,5 +469,5 @@
                 ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
             }else{
-                bedslope_work = 0.;//-pressure_flux*length ;//g*length*( 0.*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
+                //bedslope_work = 0.;//-pressure_flux*length ;//g*length*( 0.*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
                 xmom_explicit_update[k] *= 0.;
                 ymom_explicit_update[k] *= 0.;
@@ -524,5 +524,5 @@
                     ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
                 }else{
-                    bedslope_work = 0.;//-pressure_flux*length;//g*length*(0.*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
+                    //bedslope_work = 0.;//-pressure_flux*length;//g*length*(0.*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
                     xmom_explicit_update[n] *= 0.;
                     ymom_explicit_update[n] *= 0.;
@@ -633,5 +633,5 @@
 //
     free(max_bed_edgevalue);
-    free(min_bed_edgevalue);
+    //free(min_bed_edgevalue);
 
     return timestep;
@@ -671,21 +671,16 @@
 
       if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
-        
         // Set momentum to zero and ensure h is non negative
-        // NOTE: THIS IS IMPORTANT -- WE ARE SETTING MOMENTUM TO ZERO
-        //if(hc<=epsilon){
             xmomc[k] = 0.0;
             ymomc[k] = 0.0;
-        //}
 
         if (hc <= 0.0){
              // Definine the minimum bed edge value on triangle k.
-             // Setting = minimum edge value can lead to mass conservation problems
              //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
              //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
-             // Setting = minimum vertex value seems better, but might tend to be less smooth 
              //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
              bmin=zc[k]-minimum_allowed_height;
              // Minimum allowed stage = bmin
+
              // WARNING: ADDING MASS if wc[k]<bmin
              if(wc[k]<bmin){
@@ -697,6 +692,6 @@
              
 
-                 // Set vertex values as well. Seems that this shouldn't be
-                 // needed. However from memory this is important at the first
+                 // FIXME: Set vertex values as well. Seems that this shouldn't be
+                 // needed. However, from memory this is important at the first
                  // time step, for 'dry' areas where the designated stage is
                  // less than the bed centroid value
@@ -813,5 +808,5 @@
   
   double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store, *min_elevation_edgevalue, *max_elevation_edgevalue;
-  double *neighbour_wetness_scale;
+  double *cell_wetness_scale;
   int *count_wet_neighbours;
 
@@ -821,7 +816,7 @@
   ymom_centroid_store = malloc(number_of_elements*sizeof(double));
   stage_centroid_store = malloc(number_of_elements*sizeof(double));
-  min_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
+  //min_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
   max_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
-  neighbour_wetness_scale = malloc(number_of_elements*sizeof(double));
+  cell_wetness_scale = malloc(number_of_elements*sizeof(double));
   count_wet_neighbours = malloc(number_of_elements*sizeof(int));
  
@@ -838,7 +833,7 @@
           ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
 
-          min_elevation_edgevalue[k] = min(elevation_edge_values[3*k], 
-                                           min(elevation_edge_values[3*k+1],
-                                               elevation_edge_values[3*k+2]));
+          //min_elevation_edgevalue[k] = min(elevation_edge_values[3*k], 
+          //                                 min(elevation_edge_values[3*k+1],
+          //                                     elevation_edge_values[3*k+2]));
           max_elevation_edgevalue[k] = max(elevation_edge_values[3*k], 
                                            max(elevation_edge_values[3*k+1],
@@ -852,8 +847,8 @@
   for(k=0; k<number_of_elements;k++){ 
       count_wet_neighbours[k]=0;
-      neighbour_wetness_scale[k] = 0.;
-      // Cell k is wet
-      if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
-          neighbour_wetness_scale[k] = 1.;  
+      cell_wetness_scale[k] = 0.;
+      // Check if cell k is wet
+      if(stage_centroid_values[k] > elevation_centroid_values[k] + minimum_allowed_height+epsilon){
+          cell_wetness_scale[k] = 1.;  
       }
       // Count the wet neighbours
@@ -871,7 +866,5 @@
   for(k_wetdry=0; k_wetdry<2; k_wetdry++){
       //if(((stage_centroid_values[k] > max_elevation_edgevalue[k]))==k_wetdry){
-      if(((stage_centroid_values[k] < elevation_centroid_values[k]+minimum_allowed_height))==k_wetdry){
-        continue;
-      }
+      //if(((stage_centroid_values[k] < elevation_centroid_values[k]+minimum_allowed_height))==k_wetdry){
       // For partially wet cells, we now know that the edges of neighbouring
       // fully wet cells have been defined
@@ -880,4 +873,13 @@
       for (k = 0; k < number_of_elements; k++) 
       {
+
+        // First treat all 'fully wet' cells, then treat 'partially dry' cells
+        // For partially wet cells, we now know that the edges of neighbouring
+        // fully wet cells have been defined
+        if( cell_wetness_scale[k]==(1.0*k_wetdry) ){
+          continue;
+        }
+
+        // Useful indices
         k3=k*3;
         k6=k*6;
@@ -950,13 +952,16 @@
           // used
           //if(stage_centroid_values[k2]<=max_elevation_edgevalue[k2]+epsilon){
-          if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
+          //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
+          if(cell_wetness_scale[k2]==0.0){
               k2 = k ;
           }
           //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){
-          if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
+          //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
+          if(cell_wetness_scale[k0]==0.0){
               k0 = k ;
           }
           //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){
-          if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
+          //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
+          if(cell_wetness_scale[k1]==0.0){
               k1 = k ;
           }
@@ -999,86 +1004,107 @@
            
           // Treat triangles with zero or 1 wet neighbours. 
-          if ((area2 <= 0)) //|(count_wet_neighbours[k]==0))
+          if ((area2 <= 0)|(cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
           {
-              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];
-                  //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];
-
-                  // Isolated wet cell -- constant depth extrapolation, but don't let water flow uphill
-                  //stage_edge_values[k3]   = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3]);
-                  //stage_edge_values[k3+1] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1]);
-                  //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{
-                  // 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(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
-                  if(k_wetdry==0){
-                  //    // Cell is wet
-                  //    //if(count_wet_neighbours[ktmp]>0){
-                  //    //    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{
-                  //        // Constant depth extrapolation
-                      //if(stage_centroid_values[k]<stage_centroid_values[ktmp]){
-                          //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];
-                      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[ktmp];
-                      //    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]);
-                      //}
-
-                  }else{
-                      
-                  //    // This cell is 'dry'.
-                  //    // Extrapolate using the neighbouring wet cell if appropriate
-                  //    // This needs to preserve a wet-dry interface
-                      stage_edge_values[k3]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
-                      stage_edge_values[k3+1]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
-                      stage_edge_values[k3+2]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
-                  //    stage_edge_values[k3]= stage_centroid_values[ktmp] ;
-                  //    stage_edge_values[k3+1]= stage_centroid_values[ktmp];
-                  //    stage_edge_values[k3+2]= stage_centroid_values[ktmp];
-                  //    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((cell_wetness_scale[k]==1.0) | (k==k0)&&(k==k1)&&(k==k2)){
+              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{
+              // Dry cell with a single 'wet' neighbour
+
+              // Compute indices of neighbouring wet cell / edge
+              //ktmp = k0*(1-(k0==k)) + k1*(1-(k1==k)) + k2*(1-(k2==k));
+              //ii = 0*(1-(k0==k)) + 1*(1-(k1==k)) + 2*(1-(k2==k));
+              //if((ii<0)|(ii>2)){
+              //  printf("Invalid ii value \n");
+              //}
+            
+              //ktmp = 3*ktmp+neighbour_edges[3*k+ii];
+              //stage_edge_values[k3]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
+              //stage_edge_values[k3+1]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
+              //stage_edge_values[k3+2]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
+
+            //}   
+              //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];
+              //    //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];
+
+              //    // Isolated wet cell -- constant depth extrapolation, but don't let water flow uphill
+              //    //stage_edge_values[k3]   = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3]);
+              //    //stage_edge_values[k3+1] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1]);
+              //    //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{
+              //    // 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(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
+              //    if(k_wetdry==0){
+              //    //    // Cell is wet
+              //    //    //if(count_wet_neighbours[ktmp]>0){
+              //    //    //    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{
+              //    //        // Constant depth extrapolation
+              //        //if(stage_centroid_values[k]<stage_centroid_values[ktmp]){
+              //          //  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];
+              //        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[ktmp];
+              //           // 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]);
+              //        //}
+
+              //    }else{
+              //        
+              //    //    // This cell is 'dry'.
+              //    //    // Extrapolate using the neighbouring wet cell if appropriate
+              //    //    // This needs to preserve a wet-dry interface
+              //        //ktmp = 3*ktmp+neighbour_edges[3*k+ii];
+              //        //stage_edge_values[k3]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
+              //        //stage_edge_values[k3+1]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
+              //        //stage_edge_values[k3+2]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
+              //    //    stage_edge_values[k3]= stage_centroid_values[ktmp] ;
+              //    //    stage_edge_values[k3+1]= stage_centroid_values[ktmp];
+              //    //    stage_edge_values[k3+2]= stage_centroid_values[ktmp];
+              //        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];
+              //    }
+              //
+              //}
               // First order momentum / velocity extrapolation
               //stage_edge_values[k3]    = stage_centroid_values[k];
@@ -1147,5 +1173,5 @@
           // Limit the gradient
           // This is more complex for stage than other quantities
-          if((count_wet_neighbours[k]>0)){ //&&(stage_centroid_values[k] > elevation_centroid_values[k])){
+          if((count_wet_neighbours[k]>0)){//&&(k_wetdry==0)){
               limit_gradient(dqv, qmin, qmax, beta_tmp);
               stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
@@ -1153,26 +1179,27 @@
               stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
           }else{
-          //if(count_wet_neighbours[k]==0){
-              //weight=max(neighbour_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
-              //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
-              //weight=min(weight, 1.0);
-              //weight==0;
-              //weight=0.;
-              // 'Shallow flow on steep slope'
+          ////if(count_wet_neighbours[k]==0){
+          //    //weight=max(cell_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
+          //    //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
+          //    //weight=min(weight, 1.0);
+          //    //weight==0;
+          //    //weight=0.;
+          //    // 'Shallow flow on steep slope'
+          //    //limit_gradient(dqv, qmin, qmax, beta_tmp);
               stage_edge_values[k3+0] =stage_centroid_values[k] + dqv[0];
               stage_edge_values[k3+1] =stage_centroid_values[k] + dqv[1];
               stage_edge_values[k3+2] =stage_centroid_values[k] + dqv[2];
 
-              // There exists a neighbouring bed cell with elevation > minimum
-              // neighbouring stage value
-
-              // We have to be careful in this situation. Limiting the stage gradient can
-              // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
-              // Previously 'balance_deep_and_shallow' was used to deal with this issue
-              //for(i=0; i<3; i++){
-              //    stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] + 
-              //                             (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
-              //                                           +elevation_edge_values[k3+i]);
-              //}
+          //    // There exists a neighbouring bed cell with elevation > minimum
+          //    // neighbouring stage value
+
+          //    // We have to be careful in this situation. Limiting the stage gradient can
+          //    // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
+          //    // Previously 'balance_deep_and_shallow' was used to deal with this issue
+          //    //for(i=0; i<3; i++){
+          //    //    stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] + 
+          //    //                             (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
+          //    //                                           +elevation_edge_values[k3+i]);
+          //    //}
           }
            
@@ -1214,5 +1241,6 @@
           //if(stagemin>=bedmax){
           if((count_wet_neighbours[k]>0) && 
-            (stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
+             (cell_wetness_scale[k]==1.0)){
+            //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
             limit_gradient(dqv, qmin, qmax, beta_tmp);
           }else{
@@ -1266,5 +1294,6 @@
           //if(stagemin>=bedmax){
           if((count_wet_neighbours[k]>0)&&
-            (stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
+             (cell_wetness_scale[k]==1.0)){
+            //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
             limit_gradient(dqv, qmin, qmax, beta_tmp);
           }else{
@@ -1457,10 +1486,11 @@
   //    // Hack for 'dry' cells
   //    // Make sure edge value is not greater than neighbour edge value
-  //    if(stage_centroid_values[k] - elevation_centroid_values[k] < minimum_allowed_height+epsilon){
+  //    //if(stage_centroid_values[k] - elevation_centroid_values[k] < minimum_allowed_height+epsilon){
+  //    if(cell_wetness_scale[k]==0.0){
   //      for(i=0; i<3;i++){
   //          if(surrogate_neighbours[3*k+i] >= 0){
   //              ktmp=3*surrogate_neighbours[3*k+i] + neighbour_edges[3*k+i];
-  //          //if(neighbour_wetness_scale[surrogate_neighbours[3*k+i]]>0.){
-  //              stage_edge_values[3*k+i] = max(stage_edge_values[3*k+i], stage_edge_values[ktmp]);
+  //          //if(cell_wetness_scale[surrogate_neighbours[3*k+i]]>0.){
+  //              stage_edge_values[3*k+i] = min(stage_edge_values[3*k+i], stage_edge_values[ktmp]);
   //          }
   //      }
@@ -1518,7 +1548,7 @@
   free(ymom_centroid_store);
   free(stage_centroid_store);
-  free(min_elevation_edgevalue);
+  //free(min_elevation_edgevalue);
   free(max_elevation_edgevalue);
-  free(neighbour_wetness_scale);
+  free(cell_wetness_scale);
   free(count_wet_neighbours);
 
Index: /trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
===================================================================
--- /trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 8548)
+++ /trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 8549)
@@ -37,5 +37,5 @@
 
 def stagefun(x,y):
-    stge=-0.2*scale_me # +0.01*(x>0.9) 
+    stge=-0.2*scale_me #+0.01*(x>0.9) 
     #topo=topography(x,y) 
     return stge#*(stge>topo) + (topo)*(stge<=topo)
@@ -72,5 +72,5 @@
 # Associate boundary tags with boundary objects
 #----------------------------------------------
-domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
+domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
 
 #------------------------------
@@ -78,5 +78,5 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.2,finaltime=20.0):
+for t in domain.evolve(yieldstep=0.2,finaltime=40.0):
     print domain.timestepping_statistics()
     xx = domain.quantities['xmomentum'].centroid_values
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 8548)
+++ /trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py	(revision 8549)
@@ -30,14 +30,14 @@
 #------------------------------------------------------------------------------
 points, vertices, boundary = rectangular_cross(10, 10,
-len1=100.0, len2=100.0) # Mesh
+                                len1=100.0, len2=100.0) # Mesh
 domain = Domain(points, vertices, boundary) # Create domain
 domain.set_name('channel_SU_2_v2') # Output name
-domain.set_store_vertices_uniquely(True)
-domain.CFL=1.0
+#domain.set_store_vertices_uniquely(True)
+#domain.CFL=1.0
 #------------------------------------------------------------------------------
 # Setup initial conditions
 #------------------------------------------------------------------------------
 def topography(x, y):
-	return -x/10. + numpy.sin(x/10.) +abs(y-50.)/100. # linear bed slope
+	return -x/10. #+ numpy.sin(x/10.) +abs(y-50.)/100. # linear bed slope
 
 def stagetopo(x,y):
@@ -53,4 +53,5 @@
 domain.set_quantity('friction', 0.03) # Constant friction
 domain.set_quantity('stage', stagetopo)
+
 #------------------------------------------------------------------------------
 # Setup boundary conditions
@@ -69,7 +70,19 @@
 domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
 
+
+myindex=((domain.centroid_coordinates[:,0] - 50.)**2 + (domain.centroid_coordinates[:,1] - 50.)**2).argmin()
+
 #------------------------------------------------------------------------------
 # Evolve system through time
 #------------------------------------------------------------------------------
 for t in domain.evolve(yieldstep=4.0, finaltime=400.0):
-	print domain.timestepping_statistics()
+    print domain.timestepping_statistics()
+    #print domain.quantities['stage'].centroid_values[myindex] - domain.quantities['elevation'].centroid_values[myindex]
+    s3 = domain.get_flow_through_cross_section([[30., 0.0], [30., 100.]])
+    s4 = domain.get_flow_through_cross_section([[32., 0.0], [32., 100.]])
+    s5 = domain.get_flow_through_cross_section([[34., 0.0], [34., 100.]])
+    s2 = domain.get_flow_through_cross_section([[45., 0.0], [45., 100.]])
+    s1 = domain.get_flow_through_cross_section([[53., 0.0], [53., 100.]])
+    s0 = domain.get_flow_through_cross_section([[60., 0.0], [60., 100.]])
+    print 'Xsectional flow:', s0, s1, s2, s3, s4, s5
+    
