Index: trunk/anuga_work/development/gareth/experimental/bal_and/plot_utils.py
===================================================================
--- trunk/anuga_work/development/gareth/experimental/bal_and/plot_utils.py	(revision 8990)
+++ trunk/anuga_work/development/gareth/experimental/bal_and/plot_utils.py	(revision 8990)
@@ -0,0 +1,479 @@
+""" Random utilities for reading sww file data and for plotting (in ipython, or
+    in scripts)
+
+    Functionality of note:
+
+    util.get_outputs -- read the data from a single sww file into a single object
+    
+    util.combine_outputs -- read the data from a list of sww files into a single object
+    
+    util.near_transect -- for finding the indices of points 'near' to a given
+                       line, and assigning these points a coordinate along that line. 
+                       This is useful for plotting outputs which are 'almost' along a
+                       transect (e.g. a channel cross-section) -- see example below
+
+    util.sort_sww_filenames -- match sww filenames by a wildcard, and order
+                               them according to their 'time'. This means that
+                               they can be stuck together using 'combine_outputs' correctly
+
+    util.triangle_areas -- compute the areas of every triangle in a get_outputs object [ must be vertex-based]
+
+    #FIXME: util.water_volume -- compute the water volume at every time step in an sww file (needs both vertex and centroid value input). 
+ 
+    Here is an example ipython session which uses some of these functions:
+
+    > import util
+    > from matplotlib import pyplot as pyplot
+    > p=util.get_output('myfile.sww',minimum_allowed_height=0.01)
+    > p2=util.get_centroids(p,velocity_extrapolation=True)
+    > xxx=util.near_transect(p,[95., 85.], [120.,68.],tol=2.) # Could equally well use p2
+    > pyplot.ion() # Interactive plotting
+    > pyplot.scatter(xxx[1],p.vel[140,xxx[0]],color='red') # Plot along the transect
+
+    FIXME: TODO -- Convert to a single function 'get_output', which can either take a
+          single filename, a list of filenames, or a wildcard defining a number of
+          filenames, and ensure that in each case, the output will be as desired.
+
+"""
+from anuga.file.netcdf import NetCDFFile
+import numpy
+
+class combine_outputs:
+    """
+    Read in a list of filenames, and combine all their outputs into a single object.
+    e.g.:
+
+    p = util.combine_outputs(['file1.sww', 'file1_time_10000.sww', 'file1_time_20000.sww'], 0.01)
+    
+    will make an object p which has components p.x,p.y,p.time,p.stage, .... etc,
+    where the values of stage / momentum / velocity from the sww files are concatenated as appropriate.
+
+    This is nice for interactive interrogation of model outputs, or for sticking together outputs in scripts
+   
+    WARNING: It is easy to use lots of memory, if the sww files are large.
+
+    Note: If you want the centroid values, then you could subsequently use:
+
+    p2 = util.get_centroids(p,velocity_extrapolation=False)
+
+    which would make an object p2 that is like p, but holds information at centroids
+    """
+    def __init__(self, filename_list, minimum_allowed_height=1.0e-03):
+        #
+        # Go through the sww files in 'filename_list', and combine them into one object.
+        #
+
+        for i, filename in enumerate(filename_list):
+            print i, filename
+            # Store output from filename
+            p_tmp = get_output(filename, minimum_allowed_height)
+            if(i==0):
+                # Create self
+                p1=p_tmp
+            else:
+                # Append extra data to self
+                # Note that p1.x, p1.y, p1.vols, p1.elev should not change
+                assert (p1.x == p_tmp.x).all()
+                assert (p1.y == p_tmp.y).all()
+                assert (p1.vols ==p_tmp.vols).all()
+                p1.time = numpy.append(p1.time, p_tmp.time)
+                p1.stage = numpy.append(p1.stage, p_tmp.stage, axis=0)
+                p1.height = numpy.append(p1.stage, p_tmp.stage, axis=0)
+                p1.xmom = numpy.append(p1.xmom, p_tmp.xmom, axis=0)
+                p1.ymom = numpy.append(p1.ymom, p_tmp.ymom, axis=0)
+                p1.xvel = numpy.append(p1.xvel, p_tmp.xvel, axis=0)
+                p1.yvel = numpy.append(p1.yvel, p_tmp.yvel, axis=0)
+                p1.vel = numpy.append(p1.vel, p_tmp.vel, axis=0)
+
+        self.x, self.y, self.time, self.vols, self.elev, self.stage, self.xmom, \
+                self.ymom, self.xvel, self.yvel, self.vel, self.minimum_allowed_height = \
+                p1.x, p1.y, p1.time, p1.vols, p1.elev, p1.stage, p1.xmom, p1.ymom, p1.xvel,\
+                p1.yvel, p1.vel, p1.minimum_allowed_height
+
+####################
+
+def sort_sww_filenames(sww_wildcard):
+    # Function to take a 'wildcard' sww filename, 
+    # and return a list of all filenames of this type,
+    # sorted by their time.
+    # This can then be used efficiently in 'combine_outputs'
+    # if you have many filenames starting with the same pattern
+    import glob
+    filenames=glob.glob(sww_wildcard)
+    
+    # Extract time from filenames
+    file_time=range(len(filenames)) # Predefine
+     
+    for i,filename in enumerate(filenames):
+        filesplit=filename.rsplit('_time_')
+        if(len(filesplit)>1):
+            file_time[i]=int(filesplit[1].split('_0.sww')[0])
+        else:
+            file_time[i]=0         
+    
+    name_and_time=zip(file_time,filenames)
+    name_and_time.sort() # Sort by file_time
+    
+    output_times, output_names = zip(*name_and_time)
+    
+    return list(output_names)
+
+##############
+
+class get_output:
+    """Read in data from an .sww file in a convenient form
+       e.g. 
+        p = util.get_output('channel3.sww', minimum_allowed_height=0.01)
+        
+       p then contains most relevant information as e.g., p.stage, p.elev, p.xmom, etc 
+    """
+    def __init__(self, filename, minimum_allowed_height=1.0e-03):
+        self.x, self.y, self.time, self.vols, self.stage, \
+                self.height, self.elev, self.xmom, self.ymom, \
+                self.xvel, self.yvel, self.vel, self.minimum_allowed_height = \
+                read_output(filename, minimum_allowed_height)
+
+
+def read_output(filename, minimum_allowed_height):
+    # Input: The name of an .sww file to read data from,
+    #                    e.g. read_sww('channel3.sww')
+    #
+    # Purpose: To read the sww file, and output a number of variables as arrays that 
+    #          we can then manipulate (e.g. plot, interrogate)
+    #
+    # Output: x, y, time, stage, elev, xmom, ymom, xvel, yvel, vel
+
+    # Import modules
+
+
+
+    # Open ncdf connection
+    fid=NetCDFFile(filename)
+
+
+    # Read variables
+    x=fid.variables['x'][:]
+    
+    y=fid.variables['y'][:]
+
+    stage=fid.variables['stage'][:]
+
+    height=fid.variables['height'][:]
+
+    elev=fid.variables['elevation'][:]
+
+    xmom=fid.variables['xmomentum'][:]
+
+    ymom=fid.variables['ymomentum'][:]
+
+    time=fid.variables['time'][:]
+
+    vols=fid.variables['volumes'][:]
+
+
+    # Calculate velocity = momentum/depth
+    xvel=xmom*0.0
+    yvel=ymom*0.0
+
+    for i in range(xmom.shape[0]):
+        xvel[i,:]=xmom[i,:]/(height[i,:]+1.0e-06)*(height[i,:]>minimum_allowed_height)
+        yvel[i,:]=ymom[i,:]/(height[i,:]+1.0e-06)*(height[i,:]>minimum_allowed_height)
+
+    vel = (xvel**2+yvel**2)**0.5
+
+    return x, y, time, vols, stage, height, elev, xmom, ymom, xvel, yvel, vel, minimum_allowed_height
+
+##############
+
+
+class get_centroids:
+    """
+    Extract centroid values from the output of get_output.  
+    e.g.
+        p = util.get_output('my_sww.sww', minimum_allowed_height=0.01) # vertex values
+        pc=util.get_centroids(p, velocity_extrapolation=True) # centroid values
+    """
+    def __init__(self,p, velocity_extrapolation=False):
+        
+        self.time, self.x, self.y, self.stage, self.xmom,\
+             self.ymom, self.height, self.elev, self.xvel, \
+             self.yvel, self.vel= \
+             get_centroid_values(p, velocity_extrapolation)
+                                 
+
+def get_centroid_values(p, velocity_extrapolation):
+    # Input: p is the result of e.g. p=util.get_output('mysww.sww'). See the get_output class defined above
+    # Output: Values of x, y, Stage, xmom, ymom, elev, xvel, yvel, vel at centroids
+    #import numpy
+    
+    # Make 3 arrays, each containing one index of a vertex of every triangle.
+    l=len(p.vols)
+    #vols0=numpy.zeros(l, dtype='int')
+    #vols1=numpy.zeros(l, dtype='int')
+    #vols2=numpy.zeros(l, dtype='int')
+
+    # FIXME: 22/2/12/ - I think this loop is slow, should be able to do this
+    # another way
+#    for i in range(l):
+#        vols0[i]=p.vols[i][0]
+#        vols1[i]=p.vols[i][1]
+#        vols2[i]=p.vols[i][2]
+
+
+
+    vols0=p.vols[:,0]
+    vols1=p.vols[:,1]
+    vols2=p.vols[:,2]
+
+    #print vols0.shape
+    #print p.vols.shape
+
+    # Then use these to compute centroid averages 
+    x_cent=(p.x[vols0]+p.x[vols1]+p.x[vols2])/3.0
+    y_cent=(p.y[vols0]+p.y[vols1]+p.y[vols2])/3.0
+
+    stage_cent=(p.stage[:,vols0]+p.stage[:,vols1]+p.stage[:,vols2])/3.0
+    height_cent=(p.height[:,vols0]+p.height[:,vols1]+p.height[:,vols2])/3.0
+    elev_cent=(p.elev[:,vols0]+p.elev[:,vols1]+p.elev[:,vols2])/3.0
+
+    # Here, we need to treat differently the case of momentum extrapolation or
+    # velocity extrapolation
+    if velocity_extrapolation:
+        xvel_cent=(p.xvel[:,vols0]+p.xvel[:,vols1]+p.xvel[:,vols2])/3.0
+        yvel_cent=(p.yvel[:,vols0]+p.yvel[:,vols1]+p.yvel[:,vols2])/3.0
+        
+        # Now compute momenta
+        xmom_cent=stage_cent*0.0
+        ymom_cent=stage_cent*0.0
+
+        t=len(p.time)
+
+        for i in range(t):
+            xmom_cent[i,:]=xvel_cent[i,:]*(height_cent[i,:]+1e-06)*\
+                                            (height_cent[i,:]>p.minimum_allowed_height)
+            ymom_cent[i,:]=yvel_cent[i,:]*(height_cent[i,:]+1e-06)*\
+                                            (height_cent[i,:]>p.minimum_allowed_height)
+    else:
+        xmom_cent=(p.xmom[:,vols0]+p.xmom[:,vols1]+p.xmom[:,vols2])/3.0
+        ymom_cent=(p.ymom[:,vols0]+p.ymom[:,vols1]+p.ymom[:,vols2])/3.0
+
+        # Now compute velocities
+        xvel_cent=stage_cent*0.0
+        yvel_cent=stage_cent*0.0
+
+        t=len(p.time)
+
+        for i in range(t):
+            xvel_cent[i,:]=xmom_cent[i,:]/(height_cent[i,:]+1e-06)*(height_cent[i,:]>p.minimum_allowed_height)
+            yvel_cent[i,:]=ymom_cent[i,:]/(height_cent[i,:]+1e-06)*(height_cent[i,:]>p.minimum_allowed_height)
+    
+    
+    
+    
+    # Compute velocity 
+    vel_cent=(xvel_cent**2 + yvel_cent**2)**0.5
+
+    return p.time, x_cent, y_cent, stage_cent, xmom_cent,\
+             ymom_cent, height_cent, elev_cent, xvel_cent, yvel_cent, vel_cent
+
+# Make plot of stage over time.
+#pylab.close()
+#pylab.ion()
+#pylab.plot(time, stage[:,1])
+#for i in range(201):
+#    pylab.plot(time,stage[:,i])
+
+# Momentum should be 0. 
+#print 'Momentum max/min is', xmom.max() , xmom.min()
+
+#pylab.gca().set_aspect('equal')
+#pylab.scatter(x,y,c=elev,edgecolors='none')
+#pylab.colorbar()
+#
+#n=xvel.shape[0]-1
+#pylab.quiver(x,y,xvel[n,:],yvel[n,:])
+#pylab.savefig('Plot1.png')
+
+def animate_1D(time, var, x, ylab=' '): #, x=range(var.shape[1]), vmin=var.min(), vmax=var.max()):
+    # Input: time = one-dimensional time vector;
+    #        var =  array with first dimension = len(time) ;
+    #        x = (optional) vector width dimension equal to var.shape[1];
+    
+    import pylab
+    import numpy
+   
+    
+
+    pylab.close()
+    pylab.ion()
+
+    # Initial plot
+    vmin=var.min()
+    vmax=var.max()
+    line, = pylab.plot( (x.min(), x.max()), (vmin, vmax), 'o')
+
+    # Lots of plots
+    for i in range(len(time)):
+        line.set_xdata(x)
+        line.set_ydata(var[i,:])
+        pylab.draw()
+        pylab.xlabel('x')
+        pylab.ylabel(ylab)
+        pylab.title('time = ' + str(time[i]))
+    
+    return
+
+def near_transect(p, point1, point2, tol=1.):
+    # Function to get the indices of points in p less than 'tol' from the line
+    # joining (x1,y1), and (x2,y2)
+    # p comes from util.get_output('mysww.sww')
+    #
+    # e.g.
+    # import util
+    # from matplotlib import pyplot
+    # p=util.get_output('merewether_1m.sww',0.01)
+    # p2=util.get_centroids(p,velocity_extrapolation=True)
+    # #xxx=transect_interpolate.near_transect(p,[95., 85.], [120.,68.],tol=2.)
+    # xxx=util.near_transect(p,[95., 85.], [120.,68.],tol=2.)
+    # pyplot.scatter(xxx[1],p.vel[140,xxx[0]],color='red')
+
+    x1=point1[0]
+    y1=point1[1]
+
+    x2=point2[0]
+    y2=point2[1]
+
+    # Find line equation a*x + b*y + c = 0
+    # based on y=gradient*x +intercept
+    if x1!=x2:
+        gradient= (y2-y1)/(x2-x1)
+        intercept = y1 - gradient*x1
+
+        a = -gradient
+        b = 1.
+        c = -intercept
+    else:
+        #print 'FIXME: Still need to treat 0 and infinite gradients'
+        #assert 0==1
+        a=1.
+        b=0.
+        c=-x2 
+
+    # Distance formula
+    inv_denom = 1./(a**2 + b**2)**0.5
+    distp = abs(p.x*a + p.y*b + c)*inv_denom
+
+    near_points = (distp<tol).nonzero()[0]
+    
+    # Now find a 'local' coordinate for the point, projected onto the line
+    # g1 = unit vector parallel to the line
+    # g2 = vector joining (x1,y1) and (p.x,p.y)
+    g1x = x2-x1 
+    g1y = y2-y1
+    g1_norm = (g1x**2 + g1y**2)**0.5
+    g1x=g1x/g1_norm
+    g1y=g1x/g1_norm
+
+    g2x = p.x[near_points] - x1
+    g2y = p.y[near_points] - y1
+    
+    # Dot product = projected distance == a local coordinate
+    local_coord = g1x*g2x + g1y*g2y
+
+    return near_points, local_coord
+
+########################
+# TRIANGLE AREAS, WATER VOLUME
+def triangle_areas(p, subset=None):
+    # Compute areas of triangles in p -- assumes p contains vertex information
+    # subset = vector of centroid indices to include in the computation. 
+
+    if(subset is None):
+        subset=range(len(p.vols[:,0]))
+    
+    x0=p.x[p.vols[subset,0]]
+    x1=p.x[p.vols[subset,1]]
+    x2=p.x[p.vols[subset,2]]
+    
+    y0=p.y[p.vols[subset,0]]
+    y1=p.y[p.vols[subset,1]]
+    y2=p.y[p.vols[subset,2]]
+    
+    # Vectors for cross-product
+    v1_x=x0-x1
+    v1_y=y0-y1
+    #
+    v2_x=x2-x1
+    v2_y=y2-y1
+    # Area
+    area=(v1_x*v2_y-v1_y*v2_x)*0.5
+    area=abs(area)
+    return area
+
+###
+
+#def water_volume(p,p2, per_unit_area=False, subset=None):
+#    # Compute the water volume from p(vertex values) and p2(centroid values)
+#
+#    if(subset is None):
+#        subset=range(len(p2.x))
+#
+#    l=len(p2.time)
+#    area=triangle_areas(p, subset=subset)
+#    
+#    total_area=area.sum()
+#    volume=p2.time*0.
+#   
+#    # This accounts for how volume is measured in ANUGA 
+#    # Compute in 2 steps to reduce precision error (important sometimes)
+#    for i in range(l):
+#        #volume[i]=((p2.stage[i,subset]-p2.elev[subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
+#        volume[i]=((p2.stage[i,subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
+#        volume[i]=volume[i]+((-p2.elev[subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
+#    
+#    if(per_unit_area):
+#        volume=volume/total_area 
+#    
+#    return volume
+
+
+def get_triangle_containing_point(p,point):
+
+    V = p.vols
+
+    x = p.x
+    y = p.y
+
+    l = len(x)
+
+    from anuga.geometry.polygon import is_outside_polygon,is_inside_polygon
+
+    # FIXME: Horrible brute force
+    for i in xrange(l):
+        i0 = V[i,0]
+        i1 = V[i,1]
+        i2 = V[i,2]
+        poly = [ [x[i0], y[i0]], [x[i1], y[i1]], [x[i2], y[i2]] ]
+
+        if is_inside_polygon(point, poly, closed=True):
+            return i
+
+    msg = 'Point %s not found within a triangle' %str(point)
+    raise Exception(msg)
+
+
+def get_extent(p):
+
+    import numpy
+
+    x_min = numpy.min(p.x)
+    x_max = numpy.max(p.x)
+
+    y_min = numpy.min(p.y)
+    y_max = numpy.max(p.y)
+
+    return x_min, x_max, y_min, y_max
+
+
+
+
Index: trunk/anuga_work/development/gareth/experimental/bal_and/swb2_boundary_conditions.py
===================================================================
--- trunk/anuga_work/development/gareth/experimental/bal_and/swb2_boundary_conditions.py	(revision 8990)
+++ trunk/anuga_work/development/gareth/experimental/bal_and/swb2_boundary_conditions.py	(revision 8990)
@@ -0,0 +1,301 @@
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Boundary, File_boundary
+
+#from swb2_domain_ext import rotate
+#####
+
+#class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        CAN BE UNSTABLE -- CAREFUL
+#        Sets boundary values of h, uh, vh as in transmissive boundary, but 
+#        also ensures:
+#        flux[0] = 0.0,
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """Transmissive boundaries return the edge values
+#        of the volume they serve.
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#        return q
+#
+######
+#
+#class  zero_mass_flux_zero_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=0, vh=0. 
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=0.0
+#        ndamp=0.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+#class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=0, vh=vh 
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=1.0
+#        ndamp=0.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+#class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=0.0
+#        ndamp=1.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+class Characteristic_boundary(Boundary):
+    """
+
+    Example:
+
+    def waveform(t):
+        return sea_level + normalized_amplitude/cosh(t-25)**2
+
+    Bf = Nudge_boundary(domain, waveform)
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain=None, function=None):
+        """ Instantiate a
+            Nudge_boundary.
+            domain is the domain containing the boundary
+            function is the function to apply
+        """
+
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for this type boundary'
+            raise Exception, msg
+
+        if function is None:
+            msg = 'Function must be specified for this type boundary'
+            raise Exception, msg
+
+        self.domain = domain
+        self.function = function
+
+
+    def __repr__(self):
+        """ Return a representation of this instance. """
+        msg = 'Nudge_boundary'
+        msg += '(%s)' % self.domain
+        return msg
+
+
+    def evaluate(self, vol_id, edge_id):
+        """
+        """
+
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        bed = self.domain.quantities['elevation'].centroid_values[vol_id]
+        depth=max(q[0]-bed,0.0)
+        dt=self.domain.timestep
+
+        normal = self.domain.get_normal(vol_id, edge_id)
+
+
+        t = self.domain.get_time()
+
+        if hasattr(self.function, 'time'):
+            # Roll boundary over if time exceeds
+            while t > self.function.time[-1]:
+                msg = 'WARNING: domain time %.2f has exceeded' % t
+                msg += 'time provided in '
+                msg += 'Flather_boundary object.\n'
+                msg += 'I will continue, reusing the object from t==0'
+                log.critical(msg)
+                t -= self.function.time[-1]
+
+        value = self.function(t)
+        try:
+            x = float(value)
+        except:
+            x = float(value[0])
+
+        ##
+        ndotq = (normal[0]*q[1] + normal[1]*q[2]) # momentum perpendicular to the boundary
+
+        if(depth==0.):
+            q[0] = x
+            q[1] = 0. 
+            q[2] = 0. 
+
+        else:
+
+            # Asssume sub-critical flow. Set the values of the characteristics as
+            # appropriate, depending on whether we have inflow or outflow
+            tmp = (9.8/depth)**0.5
+            if(ndotq>0.):
+                # Only need to impose one characteristic
+                # Assume zero external velocity (dangerous?)
+                # FIXME: Allow external velocities to be specified in the function
+                # This would be needed to use the boundary for real inflows (e.g. run_wave.py)
+                # At the moment, it seems to work okay as a 'passive' outflow boundary (e.g. run_wave.py)
+                # Here it attenuates the wave at the outflow a bit, but doesn't do too much damage compared
+                # with the other boundary conditions
+                w1 = 0. - tmp*x
+                w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
+                w3 = ndotq/depth + tmp*q[0]
+                
+            else:
+                # Need to set 2 characteristics
+                # Assume zero external velocity (dangerous?)
+                # FIXME: Allow this to be specified with a function
+                w1 = 0. - tmp*x
+                w2 = 0.
+                #w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
+                w3 = ndotq/depth + tmp*q[0]
+
+
+            q[0] = (w3-w1)/(2*tmp)
+            qperp= (w3+w1)/2. *depth
+            qpar=  w2*depth
+
+            # So q[1], q[2] = qperp*(normal[0], normal[1]) + qpar*(-normal[1], normal[0])
+
+            q[1] = qperp*normal[0] + qpar*normal[1]
+            q[2] = qperp*normal[1] -qpar*normal[0]
+
+        return q
+###
+def rotate(q,normal,direction=1):
+    #Function for rotation -- better to mimic the shallow water one
+    # FIXME: Use C version
+    out=q*1.      
+
+    if(direction==1):
+        out[1] = q[1]*normal[0] + q[2]*normal[1]    
+        out[2] = -q[1]*normal[1] + q[2]*normal[0]
+    elif(direction==-1):
+        out[1] = q[1]*normal[0] - q[2]*normal[1]    
+        out[2] = +q[1]*normal[1] + q[2]*normal[0]
+    else:
+        print 'ERROR IN rotate'
+        assert 0==1
+
+    return out
Index: trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py
===================================================================
--- trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py	(revision 8990)
+++ trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py	(revision 8990)
@@ -0,0 +1,661 @@
+
+from anuga.shallow_water.shallow_water_domain import *
+from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain
+import numpy
+
+##############################################################################
+# Shallow Water Balanced Domain -- alternative implementation
+# 
+# FIXME: Following the methods in CITE MODSIM PAPER
+#
+##############################################################################
+
+class Domain(Sww_domain):
+
+    def __init__(self,
+                 coordinates=None,
+                 vertices=None,
+                 boundary=None,
+                 tagged_elements=None,
+                 geo_reference=None,
+                 use_inscribed_circle=False,
+                 mesh_filename=None,
+                 use_cache=False,
+                 verbose=False,
+                 full_send_dict=None,
+                 ghost_recv_dict=None,
+                 starttime=0.0,
+                 processor=0,
+                 numproc=1,
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
+
+        conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum']
+
+        evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum']         
+        other_quantities   = [ 'elevation', 'friction', 'height', 
+                               'xvelocity', 'yvelocity', 'x', 'y' ]
+
+        
+        Sww_domain.__init__(self,
+                            coordinates = coordinates,
+                            vertices = vertices,
+                            boundary = boundary,
+                            tagged_elements = tagged_elements,
+                            geo_reference = geo_reference,
+                            use_inscribed_circle = use_inscribed_circle,
+                            mesh_filename = mesh_filename,
+                            use_cache = use_cache,
+                            verbose = verbose,
+                            conserved_quantities = conserved_quantities,
+                            evolved_quantities = evolved_quantities,
+                            other_quantities = other_quantities,
+                            full_send_dict = full_send_dict,
+                            ghost_recv_dict = ghost_recv_dict,
+                            starttime = starttime,
+                            processor = processor,
+                            numproc = numproc,
+                            number_of_full_nodes = number_of_full_nodes,
+                            number_of_full_triangles = number_of_full_triangles)
+       
+        #self.forcing_terms.append(manning_friction_implicit)
+        #------------------------------------------------
+        # set some defaults
+        # Most of these override the options in config.py
+        #------------------------------------------------
+        self.set_CFL(1.0)
+        self.set_use_kinematic_viscosity(False)
+        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 
+        self.minimum_allowed_height=1.0e-03
+
+        self.use_edge_limiter=True
+        self.default_order=2
+        self.extrapolate_velocity_second_order=True
+
+        self.beta_w=1.0
+        self.beta_w_dry=0.0
+        self.beta_uh=1.0
+        self.beta_uh_dry=0.0
+        self.beta_vh=1.0
+        self.beta_vh_dry=0.0
+        
+
+        self.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 
+                 'ymomentum': 2, 'elevation': 2, 'height':2})
+
+        #self.epsilon=1.0e-100
+
+        #self.optimise_dry_cells=True
+        # "self.optimise_dry_cells=False" presently ensures that the stage is
+        # always >= minimum_bed_edge value.  Actually, we should only need to
+        # apply 'False' on the very first time step (to deal with stage values
+        # that were initialised below the bed by the user). After this, the
+        # algorithm should take care of itself, and 'True' should be okay.
+        self.optimise_dry_cells=False 
+
+        # Because gravity is treated within the flux function, 
+        # we remove it from the forcing terms.
+        #self.forcing_terms.remove(gravity)
+
+        # We need the edge_coordinates for the extrapolation
+        self.edge_coordinates=self.get_edge_midpoint_coordinates()
+
+        # We demand that vertex values are stored uniquely
+        self.set_store_vertices_smoothly(False)
+
+        self.maximum_allowed_speed=0.0
+        #self.forcing_terms.append(manning_friction_explicit)
+        #self.forcing_terms.remove(manning_friction_implicit)
+
+        # Keep track of the fluxes through the boundaries
+        self.boundary_flux_integral=numpy.ndarray(1)
+        self.boundary_flux_integral[0]=0. 
+        self.boundary_flux_sum=numpy.ndarray(1)
+        self.boundary_flux_sum[0]=0.
+
+        self.call=1 # Integer counting how many times we call compute_fluxes_central
+
+        # Integer recording the order of the time-stepping scheme
+        if(self.timestepping_method=='rk2'):
+          self.timestep_order=2
+        elif(self.timestepping_method=='euler'):
+          self.timestep_order=1
+        elif(self.timestepping_method=='rk3'):
+          self.timestep_order=3
+        else:
+          err_mess='ERROR: timestepping_method= ' + self.timestepping_method +' not supported in this solver'
+          raise Exception, err_mess
+
+        print '##########################################################################'
+        print '#'
+        print '# Using experimental audusse type solver in anuga_work/development/gareth/experimental/bal_and/'
+        print "#"
+        print "# This solver is experimental. Here are some tips on using it"
+        print "#"
+        print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values"
+        print "# , as the latter can be completely misleading near strong gradients in the flow. "
+        print "# There is a util.py script in anuga_work/development/gareth/tests which might help you extract"
+        print "# quantities at centroid values from sww files."
+        print "# Note that to accuractely compute centroid values from sww files, the files need to store "
+        print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer"
+        print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect"
+        print "# ANUGA viewer as well -- I expect this would be lots of work)"
+        print "#"
+        print '# 2) Note that many options in config.py have been overridden by the solver -- I have '
+        print '# deliberately attempted to get the solver to perform well with consistent values of '
+        print '# these parameters -- so I would advise against changing them unless you at least check that '
+        print '# it really does improve things'
+        print '#'
+        print '##########################################################################'
+
+    #-------------------------------
+    # Specialisation of Evolve
+    #-------------------------------
+    # This alters the 'evolve' method from shallow_water_domain.py so that just
+    # before the file-output step, the vertex values are replaced with the
+    # centroid values. 
+    # This is useful so that we can unambigously know what the centroid values
+    # of various parameters are
+
+    def evolve(self,
+               yieldstep=None,
+               finaltime=None,
+               duration=None,
+               skip_initial_step=False):
+        """Specialisation of basic evolve method from parent class.
+        
+            Evolve the model by 1 step.
+        """
+
+        # Call check integrity here rather than from user scripts
+        # self.check_integrity()
+
+        #print 'Into Evolve'
+
+        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
+        assert 0 <= self.beta_w <= 2.0, msg
+
+        # Initial update of vertex and edge values before any STORAGE
+        # and or visualisation.
+        # This is done again in the initialisation of the Generic_Domain
+        # evolve loop but we do it here to ensure the values are ok for storage.
+        self.distribute_to_vertices_and_edges()
+
+        if self.store is True and self.get_time() == self.get_starttime():
+            self.initialise_storage()
+        #print 'Into Generic_Domain Evolve'
+        # Call basic machinery from parent class
+        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
+                                       finaltime=finaltime, duration=duration,
+                                       skip_initial_step=skip_initial_step):
+            #print 'Out of Generic_Domain Evolve'
+            # Store model data, e.g. for subsequent visualisation
+            if self.store is True:
+                # FIXME: Extrapolation is done before writing, because I had
+                # trouble correctly computing the centroid values from the
+                # vertex values (an 'average' did not seem to work correctly in
+                # very shallow cells -- there was a discrepency between
+                # domain.quantity['blah'].centroid_values and the values
+                # computed from the sww using the vertex averge). There should
+                # be a much much more disk-efficient way to store the centroid
+                # values than this
+                self.extrapolate_second_order_edge_sw()
+
+                # Store the timestep
+                self.store_timestep()
+
+            # Pass control on to outer loop for more specific actions
+            yield(t)
+    
+    #-----------------
+    # Flux computation
+    #-----------------
+
+    ## @brief Compute fluxes and timestep suitable for all volumes in domain.
+    # @param domain The domain to calculate fluxes for.
+    def compute_fluxes(domain):
+        """Compute fluxes and timestep suitable for all volumes in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        Fluxes across each edge are scaled by edgelengths and summed up
+        Resulting flux is then scaled by area and stored in
+        explicit_update for each of the three conserved quantities
+        stage, xmomentum and ymomentum
+
+        The maximal allowable speed computed by the flux_function for each volume
+        is converted to a timestep that must not be exceeded. The minimum of
+        those is computed as the next overall timestep.
+
+        Post conditions:
+          domain.explicit_update is reset to computed flux values
+          domain.timestep is set to the largest step satisfying all volumes.
+
+        This wrapper calls the underlying C version of compute fluxes
+        """
+
+        import sys
+        from swb2_domain_ext import compute_fluxes_ext_central \
+                                      as compute_fluxes_ext
+
+        #print "."
+
+        # Shortcuts
+        #if(domain.timestepping_method!='rk2'):
+        #    err_mess='ERROR: Timestepping method is ' + domain.timestepping_method +'. '+\
+        #             'You need to use rk2 timestepping with this solver, ' +\
+        #             ' because there is presently a hack in compute fluxes central. \n'+\
+        #             ' The HACK: The timestep will only ' +\
+        #             'be recomputed on every 2nd call to compute_fluxes_central, to support'+\
+        #             ' correct treatment of wetting and drying'
+
+        #    raise Exception, err_mess 
+
+        #if(domain.timestepping_method=='rk2'):
+        #  timestep_order=2
+        #elif(domain.timestepping_method=='euler'):
+        #  timestep_order=1
+        #elif(domain.timestepping_method=='rk3'):
+        #  timestep_order=3
+        #else:
+        #  err_mess='ERROR: domain.timestepping_method= ' + domain.timestepping_method +' not supported in this solver'
+        #  raise Exception, err_mess
+
+        ##print 'timestep_order=', timestep_order
+
+        Stage = domain.quantities['stage']
+        Xmom = domain.quantities['xmomentum']
+        Ymom = domain.quantities['ymomentum']
+        Bed = domain.quantities['elevation']
+        height = domain.quantities['height']
+
+        timestep = float(sys.maxint)
+
+        domain.call+=1
+        flux_timestep = compute_fluxes_ext(timestep,
+                                           domain.epsilon,
+                                           domain.minimum_allowed_height,
+                                           domain.g,
+                                           domain.boundary_flux_sum,
+                                           domain.neighbours,
+                                           domain.neighbour_edges,
+                                           domain.normals,
+                                           domain.edgelengths,
+                                           domain.radii,
+                                           domain.areas,
+                                           domain.tri_full_flag,
+                                           Stage.edge_values,
+                                           Xmom.edge_values,
+                                           Ymom.edge_values,
+                                           Bed.edge_values,
+                                           height.edge_values,
+                                           Stage.boundary_values,
+                                           Xmom.boundary_values,
+                                           Ymom.boundary_values,
+                                           domain.boundary_flux_type,
+                                           Stage.explicit_update,
+                                           Xmom.explicit_update,
+                                           Ymom.explicit_update,
+                                           domain.already_computed_flux,
+                                           domain.max_speed,
+                                           int(domain.optimise_dry_cells),
+                                           domain.timestep_order,
+                                           Stage.centroid_values,
+                                           Xmom.centroid_values,
+                                           Ymom.centroid_values,
+                                           Bed.centroid_values, 
+                                           height.centroid_values,
+                                           Bed.vertex_values)
+
+
+        # Update the boundary flux integral
+        if(domain.timestepping_method=='rk2'):
+            if(domain.call%2==1):
+              domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
+                                                domain.boundary_flux_sum[0]*domain.timestep*0.5
+              #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
+              domain.boundary_flux_sum[0]=0.
+
+        elif(domain.timestepping_method=='euler'):
+            domain.boundary_flux_integral=0.
+            # This presently doesn't work -- this section of code may need to go elsewhere
+            #domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
+            #                                  domain.boundary_flux_sum[0]*domain.timestep
+            #domain.boundary_flux_sum[0]=0.
+
+        elif(domain.timestepping_method=='rk3'):
+            domain.boundary_flux_integral=0.
+            # This needs to be designed
+        else:
+            mess='ERROR: domain.timestepping_method', domain.timestepping_method,' method not recognised'
+            raise Exception, mess
+
+        domain.flux_timestep = flux_timestep
+
+
+
+    def protect_against_infinitesimal_and_negative_heights(domain):
+        """protect against infinitesimal heights and associated high velocities"""
+
+        from swb2_domain_ext import protect
+        #print'using swb2_protect_against ..'
+
+        # shortcuts
+        wc = domain.quantities['stage'].centroid_values
+        wv = domain.quantities['stage'].vertex_values
+        zc = domain.quantities['elevation'].centroid_values
+        zv = domain.quantities['elevation'].vertex_values
+        xmomc = domain.quantities['xmomentum'].centroid_values
+        ymomc = domain.quantities['ymomentum'].centroid_values
+        areas = domain.areas
+        xc = domain.centroid_coordinates[:,0]
+        yc = domain.centroid_coordinates[:,1] 
+
+        protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
+                domain.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas, xc, yc)
+
+    def conserved_values_to_evolved_values(self, q_cons, q_evol):
+        """Mapping between conserved quantities and the evolved quantities.
+        Used where we have a boundary condition which works with conserved
+        quantities and we now want to use them for the new well balanced
+        code using the evolved quantities
+
+        Typically the old boundary condition will set the values in q_cons,
+
+        q_evol on input will have the values of the evolved quantities at the
+        edge point (useful as it provides values for evlevation).
+
+        """
+
+        wc  = q_cons[0]
+        uhc = q_cons[1]
+        vhc = q_cons[2]
+
+
+        q_evol[0]  = wc
+        q_evol[1]  = uhc
+        q_evol[2]  = vhc
+
+        return q_evol
+
+    def distribute_to_vertices_and_edges(self):
+        """ Call correct module function """
+
+        if self.use_edge_limiter:
+            #distribute_using_edge_limiter(self)
+            self.distribute_using_edge_limiter()
+        else:
+            #distribute_using_vertex_limiter(self)
+            self.distribute_using_vertex_limiter()
+
+
+    def distribute_using_vertex_limiter(domain):
+        """Distribution from centroids to vertices specific to the SWW equation.
+
+        Precondition:
+          All quantities defined at centroids and bed elevation defined at
+          vertices.
+
+        Postcondition
+          Conserved quantities defined at vertices
+        """
+        # Remove very thin layers of water
+        domain.protect_against_infinitesimal_and_negative_heights()
+
+        # Extrapolate all conserved quantities
+        if domain.optimised_gradient_limiter:
+            # MH090605 if second order,
+            # perform the extrapolation and limiting on
+            # all of the conserved quantities
+
+            if (domain._order_ == 1):
+                for name in domain.conserved_quantities:
+                    Q = domain.quantities[name]
+                    Q.extrapolate_first_order()
+            elif domain._order_ == 2:
+                domain.extrapolate_second_order_sw()
+            else:
+                raise Exception('Unknown order')
+        else:
+            # Old code:
+            for name in domain.conserved_quantities:
+                Q = domain.quantities[name]
+
+                if domain._order_ == 1:
+                    Q.extrapolate_first_order()
+                elif domain._order_ == 2:
+                    Q.extrapolate_second_order_and_limit_by_vertex()
+                else:
+                    raise Exception('Unknown order')
+
+        # Take bed elevation into account when water heights are small
+        #balance_deep_and_shallow(domain)
+
+        # Compute edge values by interpolation
+        for name in domain.conserved_quantities:
+            Q = domain.quantities[name]
+            Q.interpolate_from_vertices_to_edges()
+
+
+    def distribute_using_edge_limiter(domain):
+        """Distribution from centroids to edges specific to the SWW eqn.
+
+        Precondition:
+          All quantities defined at centroids and bed elevation defined at
+          vertices.
+
+        Postcondition
+          Conserved quantities defined at vertices
+        """
+        #from swb2_domain import protect_against_infinitesimal_and_negative_heights
+
+        # Remove very thin layers of water
+        #print 'Before protect'
+        domain.protect_against_infinitesimal_and_negative_heights()
+        #print 'After protect'
+
+        #for name in domain.conserved_quantities:
+        #    Q = domain.quantities[name]
+        #    if domain._order_ == 1:
+        #        Q.extrapolate_first_order()
+        #    elif domain._order_ == 2:
+        #        #Q.extrapolate_second_order_and_limit_by_edge()
+        #        # FIXME: This use of 'break' is hacky
+        #        domain.extrapolate_second_order_edge_sw()
+        #        break
+        #    else:
+        #        raise Exception('Unknown order')
+        
+        #print 'Before extrapolate'
+        domain.extrapolate_second_order_edge_sw()
+        #print 'After extrapolate'
+
+        #balance_deep_and_shallow(domain)
+
+        # Compute vertex values by interpolation
+        #for name in domain.conserved_quantities:
+        #    Q = domain.quantities[name]
+        #    Q.interpolate_from_edges_to_vertices()
+        #    Q.interpolate_from_vertices_to_edges()
+
+
+    def update_other_quantities(self):
+        """ There may be a need to calculates some of the other quantities
+        based on the new values of conserved quantities
+
+        But that is not needed in this version of the solver.
+        """
+        pass
+        # The centroid values of height and x and y velocity
+        # might not have been setup
+        
+        #self.update_centroids_of_velocities_and_height()
+        #
+        #for name in ['height', 'xvelocity', 'yvelocity']:
+        #    Q = self.quantities[name]
+        #    if self._order_ == 1:
+        #        Q.extrapolate_first_order()
+        #    elif self._order_ == 2:
+        #        Q.extrapolate_second_order_and_limit_by_edge()
+        #    else:
+        #        raise Exception('Unknown order')
+
+
+    def update_centroids_of_velocities_and_height(self):
+        """Calculate the centroid values of velocities and height based
+        on the values of the quantities stage and x and y momentum
+
+        Assumes that stage and momentum are up to date
+
+        Useful for kinematic viscosity calculations
+        """
+        pass
+        ## For shallow water we need to update height xvelocity and yvelocity
+
+        ##Shortcuts
+        #W  = self.quantities['stage']
+        #UH = self.quantities['xmomentum']
+        #VH = self.quantities['ymomentum']
+        #H  = self.quantities['height']
+        #Z  = self.quantities['elevation']
+        #U  = self.quantities['xvelocity']
+        #V  = self.quantities['yvelocity']
+
+        ##print num.min(W.centroid_values)
+
+        ## Make sure boundary values of conserved quantites
+        ## are consistent with value of functions at centroids
+        ##self.distribute_to_vertices_and_edges()
+        #Z.set_boundary_values_from_edges()
+
+        ##W.set_boundary_values_from_edges()
+        ##UH.set_boundary_values_from_edges()
+        ##VH.set_boundary_values_from_edges()
+
+        ## Update height values
+        #H.set_values(W.centroid_values-Z.centroid_values, location='centroids')
+        #H.set_boundary_values( num.where(W.boundary_values-Z.boundary_values>=0,
+        #                                 W.boundary_values-Z.boundary_values, 0.0))
+
+        ##assert num.min(H.centroid_values) >= 0
+        ##assert num.min(H.boundary_values) >= 0
+
+        ##Aliases
+        #uh_C  = UH.centroid_values
+        #vh_C  = VH.centroid_values
+        #h_C   = H.centroid_values
+
+        #uh_B  = UH.boundary_values
+        #vh_B  = VH.boundary_values
+        #h_B   = H.boundary_values
+
+        #H0 = 1.0e-8
+        #
+        #U.set_values(uh_C/(h_C + H0/h_C), location='centroids')
+        #V.set_values(vh_C/(h_C + H0/h_C), location='centroids')
+
+        #U.set_boundary_values(uh_B/(h_B + H0/h_B))
+        #V.set_boundary_values(vh_B/(h_B + H0/h_B))
+
+
+
+    def extrapolate_second_order_edge_sw(self):
+        """Wrapper calling C version of extrapolate_second_order_sw, using
+            edges instead of vertices in the extrapolation step. The routine
+            then interpolates from edges to vertices. 
+            The momentum extrapolation / interpolation is based on either 
+            momentum or velocity, depending on the choice of 
+            extrapolate_velocity_second_order.
+
+        self  the domain to operate on
+
+        """
+
+        from swb2_domain_ext import extrapolate_second_order_edge_sw as extrapol2
+
+        # Shortcuts
+        Stage = self.quantities['stage']
+        Xmom = self.quantities['xmomentum']
+        Ymom = self.quantities['ymomentum']
+        Elevation = self.quantities['elevation']
+        height=self.quantities['height']
+
+        extrapol2(self,
+                  self.surrogate_neighbours,
+                  self.neighbour_edges,
+                  self.number_of_boundaries,
+                  self.centroid_coordinates,
+                  Stage.centroid_values,
+                  Xmom.centroid_values,
+                  Ymom.centroid_values,
+                  Elevation.centroid_values,
+                  height.centroid_values,
+                  self.edge_coordinates,
+                  Stage.edge_values,
+                  Xmom.edge_values,
+                  Ymom.edge_values,
+                  Elevation.edge_values,
+                  height.edge_values,
+                  Stage.vertex_values,
+                  Xmom.vertex_values,
+                  Ymom.vertex_values,
+                  Elevation.vertex_values,
+                  height.vertex_values,
+                  int(self.optimise_dry_cells),
+                  int(self.extrapolate_velocity_second_order))
+
+#    def set_boundary(self, boundary_map):
+#        ## Specialisation of set_boundary, which also updates the 'boundary_flux_type'
+#        Sww_domain.set_boundary(self,boundary_map)
+#        
+#        # Add a flag which can be used to distinguish flux boundaries within
+#        # compute_fluxes_central
+#        # Initialise to zero (which means 'not a flux_boundary')
+#        self.boundary_flux_type = self.boundary_edges*0
+#       
+#        # HACK to set the values of domain.boundary_flux 
+#        for i in range(len(self.boundary_objects)):
+#            # Record the first 10 characters of the name of the boundary.
+#            # FIXME: There must be a better way than this!
+#            bndry_name = self.boundary_objects[i][1].__repr__()[0:10]
+#
+#            if(bndry_name=='zero_mass_'):
+#                # Create flag 'boundary_flux_type' that can be read in
+#                # compute_fluxes_central
+#                self.boundary_flux_type[i]=1
+#
+#
+#def manning_friction_implicit(domain):
+#    """Apply (Manning) friction to water momentum
+#    Wrapper for c version
+#    """
+#
+#    from shallow_water_ext import manning_friction_flat
+#    from shallow_water_ext import manning_friction_sloped
+#
+#    xmom = domain.quantities['xmomentum']
+#    ymom = domain.quantities['ymomentum']
+#
+#    x = domain.get_vertex_coordinates()
+#    
+#    w = domain.quantities['stage'].centroid_values
+#    z = domain.quantities['elevation'].vertex_values
+#
+#    uh = xmom.centroid_values
+#    vh = ymom.centroid_values
+#    eta = domain.quantities['friction'].centroid_values
+#
+#    xmom_update = xmom.semi_implicit_update
+#    ymom_update = ymom.semi_implicit_update
+#
+#    eps = domain.epsilon
+#    g = domain.g
+#
+#    if domain.use_sloped_mannings:
+#        manning_friction_sloped(g, eps, x, w, uh, vh, z, eta, xmom_update, \
+#                                ymom_update)
+#    else:
+#        manning_friction_flat(g, eps, w, uh, vh, z, eta, xmom_update, \
+#                                #ymom_update)
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 8990)
+++ trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c	(revision 8990)
@@ -0,0 +1,2277 @@
+// Python - C extension module for shallow_water.py
+//
+// To compile (Python2.6):
+//  gcc -c swb2_domain_ext.c -I/usr/include/python2.6 -o domain_ext.o -Wall -O
+//  gcc -shared swb2_domain_ext.o  -o swb2_domain_ext.so
+//
+// or use python compile.py
+//
+// See the module swb_domain.py for more documentation on 
+// how to use this module
+//
+//
+// Stephen Roberts, ANU 2009
+// Ole Nielsen, GA 2004
+// Gareth Davies, GA 2011
+
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "math.h"
+#include <stdio.h>
+//#include "numpy_shim.h"
+
+// Shared code snippets
+#include "util_ext.h"
+
+
+const double pi = 3.14159265358979;
+
+// Computational function for rotation
+int _rotate(double *q, double n1, double n2) {
+  /*Rotate the last  2 coordinates of q (q[1], q[2])
+    from x,y coordinates to coordinates based on normal vector (n1, n2).
+
+    Result is returned in array 2x1 r
+    To rotate in opposite direction, call rotate with (q, n1, -n2)
+
+    Contents of q are changed by this function */
+
+
+  double q1, q2;
+
+  // Shorthands
+  q1 = q[1];  // x coordinate
+  q2 = q[2];  // y coordinate
+
+  // Rotate
+  q[1] =  n1*q1 + n2*q2;
+  q[2] = -n2*q1 + n1*q2;
+
+  return 0;
+}
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
+//static inline double _compute_speed(double *uh, 
+double _compute_speed(double *uh, 
+		      double *h, 
+		      double epsilon, 
+		      double h0,
+		      double limiting_threshold) {
+  
+  double u;
+
+  if (*h < limiting_threshold) {   
+    // Apply limiting of speeds according to the ANUGA manual
+    if (*h < epsilon) {
+      *h = max(0.0,*h);  // Could have been negative
+      //*h = 0.0;  // Could have been negative
+      u = 0.0;
+    } else {
+      u = *uh/(*h + h0/ *h);    
+    }
+  
+
+    // Adjust momentum to be consistent with speed
+    *uh = u * *h;
+  } else {
+    // We are in deep water - no need for limiting
+    u = *uh/ *h;
+  }
+  
+  return u;
+}
+
+// minmod limiter
+int _minmod(double a, double b){
+    // Compute minmod
+
+    if(sign(a)!=sign(b)){
+        return 0.0;
+    }else{
+        return min(fabs(a), fabs(b))*sign(a);
+    }
+
+
+}
+
+// Innermost flux function (using stage w=z+h)
+int _flux_function_central(double *q_left, double *q_right,
+                           double h_left, double h_right,
+                           double hle, double hre,
+                           double n1, double n2,
+                           double epsilon, 
+                           double h0,
+                           double limiting_threshold, 
+                           double g,
+                           double *edgeflux, double *max_speed, 
+                           double *pressure_flux, double hc, 
+                           double hc_n, double speed_max_last) 
+{
+
+  /*Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of the 'stage', w = h+z using
+    the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+  */
+
+  int i;
+
+  double w_left,  uh_left, vh_left, u_left;
+  double w_right, uh_right, vh_right, u_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double denom, inverse_denominator;
+  double uint, t1, t2, t3, min_speed, tmp;
+  // Workspace (allocate once, use many)
+  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
+
+  // Copy conserved quantities to protect from modification
+  q_left_rotated[0] = q_left[0];
+  q_right_rotated[0] = q_right[0];
+  q_left_rotated[1] = q_left[1];
+  q_right_rotated[1] = q_right[1];
+  q_left_rotated[2] = q_left[2];
+  q_right_rotated[2] = q_right[2];
+
+  // Align x- and y-momentum with x-axis
+  _rotate(q_left_rotated, n1, n2);
+  _rotate(q_right_rotated, n1, n2);
+
+
+  //z=max(zl,zr);
+  //z = 0.5*(z_left + z_right); // Average elevation values. 
+                              // Even though this will nominally allow 
+			      // for discontinuities in the elevation data, 
+			      // there is currently no numerical support for
+			      // this so results may be strange near 
+			      // jumps in the bed.
+
+  // Compute speeds in x-direction
+  w_left = q_left_rotated[0];          
+  //h_left = max(w_left - z,0.);
+  //h_left= max(hle+zl-z,0.);
+  uh_left = q_left_rotated[1];
+  if(h_left>h0){
+    u_left = uh_left/max(h_left, 1.0e-06);
+    uh_left=h_left*u_left;
+  }else{
+    u_left=0.;
+    uh_left=h_left*u_left;
+  }
+  //u_left = _compute_speed(&uh_left, &h_left, 
+  //			  epsilon, h0, limiting_threshold);
+
+  w_right = q_right_rotated[0];
+
+  //h_right = max(w_right - z, 0.);
+  //h_right= max(hre+zr-z,0.);
+
+  uh_right = q_right_rotated[1];
+  //u_right = _compute_speed(&uh_right, &h_right, 
+  //			   epsilon, h0, limiting_threshold);
+  if(h_right>h0){
+    u_right = uh_right/max(h_right, 1.0e-06);
+    uh_right=h_right*u_right;
+  }else{
+    u_right=0.;
+    uh_right=h_right*u_right;
+  }
+
+  // Momentum in y-direction
+  vh_left  = q_left_rotated[2];
+  vh_right = q_right_rotated[2];
+
+  // Limit y-momentum if necessary 
+  // Leaving this out, improves speed significantly (Ole 27/5/2009)
+  // All validation tests pass, so do we really need it anymore? 
+  //_compute_speed(&vh_left, &h_left, 
+  //		 epsilon, h0, limiting_threshold);
+  //_compute_speed(&vh_right, &h_right, 
+  //		 epsilon, h0, limiting_threshold);
+
+  // Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left);
+  soundspeed_right = sqrt(g*h_right);  
+  
+  // Code to use fast square root optimisation if desired.
+  // Timings on AMD 64 for the Okushiri profile gave the following timings
+  //
+  // SQRT           Total    Flux
+  //=============================
+  //
+  // Ref            405s     152s
+  // Fast (dbl)     453s     173s
+  // Fast (sng)     437s     171s
+  //
+  // Consequently, there is currently (14/5/2009) no reason to use this 
+  // approximation.
+  
+  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
+  //soundspeed_right = fast_squareroot_approximation(g*h_right);
+
+  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
+  if (s_max < 0.0) 
+  {
+    s_max = 0.0;
+  }
+
+  if( hc < h0){
+    s_max = 0.0;
+  }
+
+
+  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
+  if (s_min > 0.0)
+  {
+    s_min = 0.0;
+  }
+  
+  if( hc_n < h0){
+    s_min = 0.0;
+  }
+  
+  // Flux formulas
+  flux_left[0] = u_left*h_left;
+  //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;
+  //}else{
+  //  flux_left[1] = 0.;//u_left*uh_left; //+ 0.5*g*h_left*h_left;
+  //  flux_left[2] = 0.;//u_left*vh_left;
+  //}  
+
+  flux_right[0] = u_right*h_right;
+  //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;
+  //}else{
+  //  flux_right[1] = 0.; //u_right*uh_right ; //+ 0.5*g*h_right*h_right;
+  //  flux_right[2] = 0.; //u_right*vh_right;
+  //}
+
+  // Flux computation
+  denom = s_max - s_min;
+  //if (denom < epsilon) 
+  if (0==1) 
+  { 
+    // Both wave speeds are very small
+    memset(edgeflux, 0, 3*sizeof(double));
+
+    *max_speed = 0.0;
+    //*pressure_flux = 0.0;
+    *pressure_flux = 0.5*g*0.5*(h_left*h_left+h_right*h_right);
+  } 
+  else 
+  {
+    // Maximal wavespeed
+    *max_speed = max(s_max, -s_min);
+    //min_speed=min(s_max, -s_min);
+  
+    inverse_denominator = 1.0/max(denom,1.0e-100);
+    for (i = 0; i < 3; i++) 
+    {
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+
+      //if(i==0){
+      //    // Standard smoothing term
+      //    edgeflux[i] += (s_max*s_min)*(h_right - h_left);
+      //}else{
+          // Standard smoothing term
+      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
+      //}
+
+      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;
+    
+
+    // Rotate back
+    _rotate(edgeflux, n1, -n2);
+  }
+  
+  return 0;
+}
+
+
+// Computational function for flux computation
+double _compute_fluxes_central(int number_of_elements,
+        double timestep,
+        double epsilon,
+        double H0,
+        double g,
+        double* boundary_flux_sum,
+        long* neighbours,
+        long* neighbour_edges,
+        double* normals,
+        double* edgelengths,
+        double* radii,
+        double* areas,
+        long* tri_full_flag,
+        double* stage_edge_values,
+        double* xmom_edge_values,
+        double* ymom_edge_values,
+        double* bed_edge_values,
+        double* height_edge_values,
+        double* stage_boundary_values,
+        double* xmom_boundary_values,
+        double* ymom_boundary_values,
+        long*   boundary_flux_type,
+        double* stage_explicit_update,
+        double* xmom_explicit_update,
+        double* ymom_explicit_update,
+        long* already_computed_flux,
+        double* max_speed_array,
+        int optimise_dry_cells, 
+        int timestep_order,
+        double* stage_centroid_values,
+        double* xmom_centroid_values,
+        double* ymom_centroid_values,
+        double* bed_centroid_values,
+        double* height_centroid_values,
+        double* bed_vertex_values) {
+    // Local variables
+    double max_speed, length, inv_area, zl, zr;
+    double h0 = H0*2.0;//H0*H0;//H0*H0; // This ensures a good balance when h approaches H0.
+    double h_left, h_right, z_half ;  // For andusse scheme
+
+    double limiting_threshold = 10 * H0 ;//H0; //10 * H0; // Avoid applying limiter below this
+    // threshold for performance reasons.
+    // See ANUGA manual under flux limiting
+    int k, i, m, n,j, ii;
+    int ki, nm = 0, ki2,ki3, nm3; // Index shorthands
+
+    // Workspace (making them static actually made function slightly slower (Ole))
+    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
+    double stage_edges[3];//Work array
+    double bedslope_work, local_timestep;
+    int neighbours_wet[3];//Work array
+    int useint;
+    double hle, hre, zc, zc_n; 
+    double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
+    static long call = 1; // Static local variable flagging already computed flux
+    double speed_max_last, vol, smooth;
+
+    double *count_wet_edges, *edgeflux_store, *pressuregrad_store;
+
+    //count_wet_neighbours = malloc(number_of_elements*sizeof(double));
+    count_wet_edges = malloc(number_of_elements*sizeof(double));
+    edgeflux_store = malloc(number_of_elements*9*sizeof(double));
+    pressuregrad_store = malloc(number_of_elements*3*sizeof(double));
+
+    call++; // Flag 'id' of flux calculation for this timestep
+
+    // Set explicit_update to zero for all conserved_quantities.
+    // This assumes compute_fluxes called before forcing terms
+    memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double));
+    memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
+    memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
+
+    local_timestep=timestep;
+
+    //printf("timestep_order %lf \n", timestep_order*1.0);
+
+    // Compute previous call max_speed
+    speed_max_last=0.0;
+    for (k = 0; k < number_of_elements; k++){
+        
+        speed_max_last=max(speed_max_last, max_speed_array[k]);
+    }
+
+    // For all triangles
+    for (k = 0; k < number_of_elements; k++) {
+        // Loop through neighbours and compute edge flux for each
+        for (i = 0; i < 3; i++) {
+            ki = k * 3 + i; // Linear index to edge i of triangle k
+            ki2 = 2 * ki; //k*6 + i*2
+            ki3 = 3*ki; 
+
+            if (already_computed_flux[ki] == call) {
+                // We've already computed the flux across this edge
+                continue;
+            }
+
+            // Get left hand side values from triangle k, edge i
+            ql[0] = stage_edge_values[ki];
+            ql[1] = xmom_edge_values[ki];
+            ql[2] = ymom_edge_values[ki];
+            zl = bed_edge_values[ki];
+            hc = height_centroid_values[k];
+            zc = bed_centroid_values[k];
+            hle= height_edge_values[ki];
+
+            // Get right hand side values either from neighbouring triangle
+            // or from boundary array (Quantities at neighbour on nearest face).
+            n = neighbours[ki];
+            hc_n = hc;
+            zc_n = bed_centroid_values[k];
+            if (n < 0) {
+                // Neighbour is a boundary condition
+                m = -n - 1; // Convert negative flag to boundary index
+
+                qr[0] = stage_boundary_values[m];
+                qr[1] = xmom_boundary_values[m];
+                qr[2] = ymom_boundary_values[m];
+                zr = zl; // Extend bed elevation to boundary
+                hre= hle; //-qr[0]+ql[0];
+                //hc_n = hc_n + qr[0] - ql[0];
+                //hc_n = qr[0] - zr;
+            } else {
+                // Neighbour is a real triangle
+                hc_n = max(stage_centroid_values[n] - bed_centroid_values[n],0.0);
+                zc_n = bed_centroid_values[n];
+                m = neighbour_edges[ki];
+                nm = n * 3 + m; // Linear index (triangle n, edge m)
+                nm3 = nm*3;
+
+                qr[0] = stage_edge_values[nm];
+                qr[1] = xmom_edge_values[nm];
+                qr[2] = ymom_edge_values[nm];
+                zr = bed_edge_values[nm];
+                hre = height_edge_values[nm];
+                //hc_n = height_edge_values[nm]
+            }
+          
+            // Audusse magic 
+            z_half=max(zl,zr);
+            h_left= max(hle+zl-z_half,0.);
+            h_right= max(hre+zr-z_half,0.);
+            //if (fabs(zl-zr)>1.0e-10) {
+            //    report_python_error(AT,"Discontinuous Elevation");
+            //    return 0.0;
+            //}
+           
+
+            //if((count_wet_edges[k]<=1.0 | count_wet_neighbours[k]<=1.0) & 
+            //    (n>=0 && (count_wet_edges[n]<=1.0 | count_wet_neighbours[n]<=1.0))){
+            //    smooth=0.0;
+            //}
+            
+            // Edge flux computation (triangle k, edge i)
+            _flux_function_central(ql, qr, 
+                    h_left, h_right,
+                    hle, hre,
+                    normals[ki2], normals[ki2 + 1],
+                    epsilon, h0, limiting_threshold, g,
+                    edgeflux, &max_speed, &pressure_flux, hc, hc_n, 
+                    speed_max_last);
+
+            // Multiply edgeflux by edgelength
+            length = edgelengths[ki];
+            edgeflux[0] *= length;
+            edgeflux[1] *= length;
+            edgeflux[2] *= length;
+
+            edgeflux_store[ki3 + 0 ] = -edgeflux[0];
+            edgeflux_store[ki3 + 1 ] = -edgeflux[1];
+            edgeflux_store[ki3 + 2 ] = -edgeflux[2];
+
+            // Update triangle k with flux from edge i
+            // bedslope_work contains all gravity related terms -- weighting of
+            // conservative and non-conservative versions.
+
+            //bedslope_work=length*(g*0.5*(-hle*hle+h_left*h_left +(h_left+hc)*(z_half-bed_centroid_values[k]))+pressure_flux);
+            //bedslope_work=length*(g*0.5*(-hc*hc + h_left*h_left)+pressure_flux);
+            //Audusse approach
+            //bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle)+pressure_flux);
+            //bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -hle*zl)+pressure_flux);
+            bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux);
+
+            //bedslope_work=length*(g*(-0.5*max(stage_centroid_values[k] - z_half, 0.)*(stage_centroid_values[k]-z_half) + pressure_flux));
+            pressuregrad_store[ki]=bedslope_work;
+            
+            //if((k==2 | n==2) | (k==4|n==4)){
+            //if(n<0){
+            //    printf("k n: %d,%d \n", k, n);
+            //    printf("z's: %e,%e \n", zl, zr);
+            //    printf("h_LR's: %e,%e \n", h_left,h_right);
+            //    printf("h's: %e,%e \n", hle,hre);
+            //    printf("pf: %e \n", pressure_flux);
+            //    printf("flux: %e, %e, %e\n", edgeflux[0], edgeflux[1], edgeflux[2]);
+            //    printf("bs: %e \n", bedslope_work);
+            //    printf("\n");
+            //} 
+
+            ////if(hc> h0){
+            ////if(stage_centroid_values[k] > count_wet_edges[k]){
+            //  bedslope_work = length*(g*(hc*ql[0]-0.5*hle*hle) + 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)) 
+            //  //                        + pressure_flux);
+            //  //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
+            //  pressuregrad_store[ki] = bedslope_work;
+            //}else{
+            //  // 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]-0.0*hle*hle) + 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 
+            //  pressuregrad_store[ki] = bedslope_work;
+            //}
+
+            //if(count_wet_edges[k]<=1.0) pressuregrad_store[ki]=0.0;
+
+            already_computed_flux[ki] = call; // #k Done
+
+
+            // Update neighbour n with same flux but reversed sign
+            if (n >= 0) {
+
+                edgeflux_store[nm3 + 0 ] = edgeflux[0];
+                edgeflux_store[nm3 + 1 ] = edgeflux[1];
+                edgeflux_store[nm3 + 2 ] = edgeflux[2];
+
+                //bedslope_work=length*(-g*0.5*(h_right*h_right-hre*hre)+pressure_flux);
+                //bedslope_work=length*(-g*0.5*(h_right*h_right-hre*hre-hre*zr)+pressure_flux);
+                bedslope_work=length*(-g*0.5*(h_right*h_right-hre*hre-(hre+hc_n)*(zr-zc_n))+pressure_flux);
+                //bedslope_work=length*(g*(-0.5*max(stage_centroid_values[n] - z_half, 0.)*(stage_centroid_values[n]-z_half) + pressure_flux));
+                pressuregrad_store[nm]=bedslope_work;
+                //if(hc_n> h0){
+                ////if(stage_centroid_values[n] > count_wet_edges[n]){
+                //  bedslope_work = length*(g*(hc_n*qr[0]-0.5*hre*hre) + 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);
+                //  //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
+                //  pressuregrad_store[nm] = 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_n*qr[0]-0.0*hre*hre) + 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];  
+                //  pressuregrad_store[nm] = bedslope_work;
+                //}
+            
+                //if(count_wet_edges[n]<=1.0) pressuregrad_store[nm]=0.0;
+
+                already_computed_flux[nm] = call; // #n Done
+            }
+
+            //if(n==576 | n==579) pressuregrad_store[nm]=0.;
+            // Update timestep based on edge i and possibly neighbour n
+            // GD HACK:
+            // FIXME: Presently this 'hacked' to only update the timestep on
+            // every 2nd call.  This works for rk2 timestepping only, and is
+            // needed for correct wetting and drying treatment
+            if ((tri_full_flag[k] == 1) ) {
+
+                if((call-2)%timestep_order!=0) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep 
+
+                if (max_speed > epsilon) {
+                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
+
+                    // CFL for triangle k
+                    local_timestep = min(local_timestep, radii[k] / max_speed);
+
+                    if (n >= 0) {
+                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
+                        local_timestep = min(local_timestep, radii[n] / max_speed);
+                    }
+
+                }
+            }
+
+        } // End edge i (and neighbour n)
+
+        // Keep track of maximal speeds
+        max_speed_array[k] = max_speed;
+
+    } // End triangle k
+ 
+    // GD HACK 
+    // Limit edgefluxes, for mass conservation near wet/dry cells
+    for(k=0; k< number_of_elements; k++){
+            //continue;
+            hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
+            // Loop over every edge
+            for(i = 0; i<3; i++){
+                if(i==0){
+                    // Add up the outgoing flux through the cell
+                    outgoing_mass_edges=0.0;
+                    for(useint=0; useint<3; useint++){
+                        if(edgeflux_store[3*(3*k+useint)]< 0.){
+                            //outgoing_mass_edges+=1.0;
+                            outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]);
+                        }
+                    }
+                    outgoing_mass_edges*=local_timestep;
+                }
+
+                ki=3*k+i;   
+                ki2=ki*2;
+                ki3 = ki*3;
+                
+                // Prevent outflow from 'seriously' dry cells
+                // Idea: The cell will not go dry if:
+                // total_outgoing_flux <= cell volume = Area_triangle*hc
+                vol=areas[k]*hc;
+                if((edgeflux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){
+                    
+                    // This bound could be improved (e.g. we could actually sum
+                    // the + and - fluxes and check if they are too large).
+                    // However, the advantage is that we don't have to worry
+                    // about subsequent changes to the + edgeflux caused by
+                    // constraints associated with neighbouring triangles.
+                    tmp = vol/(-(outgoing_mass_edges)) ;
+                    if(tmp< 1.0){
+                        edgeflux_store[ki3+0]*=tmp;
+                        edgeflux_store[ki3+1]*=tmp;
+                        edgeflux_store[ki3+2]*=tmp;
+
+                        // Compute neighbour edge index
+                        n = neighbours[ki];
+                        if(n>=0){
+                            nm = 3*n + neighbour_edges[ki];
+                            nm3 = nm*3;
+                            edgeflux_store[nm3+0]*=tmp;
+                            edgeflux_store[nm3+1]*=tmp;
+                            edgeflux_store[nm3+2]*=tmp;
+                        }
+                    }
+                }
+        }
+     }
+
+    // Now add up stage, xmom, ymom explicit updates
+    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.;
+
+        for(i=0;i<3;i++){
+            ki=3*k+i;   
+            ki2=ki*2;
+            ki3 = ki*3;
+            n=neighbours[ki];
+
+            // GD HACK
+            // Option to limit advective fluxes
+            if(hc > -h0){
+                stage_explicit_update[k] += edgeflux_store[ki3+0];
+                // Cut these terms for shallow flows, and protect stationary states from round-off error
+                xmom_explicit_update[k] += edgeflux_store[ki3+1];
+                ymom_explicit_update[k] += edgeflux_store[ki3+2];
+            }else{
+                stage_explicit_update[k] += edgeflux_store[ki3+0];
+            }
+
+
+            // If neighbour is a boundary condition, add the flux to the boundary_flux_integral
+            if(n<0){ 
+                boundary_flux_sum[0] += edgeflux_store[ki3];
+            }
+
+            // GD HACK
+            // Compute bed slope term
+            if(hc > -h0){
+                xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki];
+                ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki];
+            }else{
+                xmom_explicit_update[k] *= 0.;
+                ymom_explicit_update[k] *= 0.;
+            }
+
+        } // end edge i
+
+        // Normalise triangle k by area and store for when all conserved
+        // quantities get updated
+        inv_area = 1.0 / areas[k];
+        stage_explicit_update[k] *= inv_area;
+        xmom_explicit_update[k] *= inv_area;
+        ymom_explicit_update[k] *= inv_area;
+    }  // end cell k
+
+    if((call-2)%timestep_order==0) timestep=local_timestep; // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step
+
+    // Look for 'dry' edges, and set momentum (& component of momentum_update)
+    // in that direction to zero. If we don't do this, then it is possible for
+    // the bedslope term normal to that edge to be nonzero, and for momentum to
+    // keep 'building' up in the cell without being able to flow out.
+    //
+    //for(k=0; k<number_of_elements; k++){
+    //    // Ignore "Dry" centroid cells, which will automatically have their
+    //    // momentum set to zero (in _protect)
+    //    if(stage_centroid_values[k] - bed_centroid_values[k] > H0){
+    //        for(i=0;i<3;i++){
+    //            ki=3*k+i;
+    //            if(stage_edge_values[ki] - bed_edge_values[ki] <=0.){
+    //              // Compute magnitude of momentum_update in direction normal to edge
+    //              tmp =  xmom_explicit_update[k]*normals[2*ki] + ymom_explicit_update[k]*normals[2*ki+1];
+    //              if(tmp > 0.){
+    //                  // Set component of momentum_update normal to edge to zero
+    //                  // [equivalently, subtract component of momentum_update normal to edge]
+    //                  xmom_explicit_update[k] -= tmp*normals[2*ki]; 
+    //                  ymom_explicit_update[k] -= tmp*normals[2*ki+1];  
+    //                  
+    //              }
+    //              
+    //              // Compute magnitude of momentum in direction normal to edge
+    //              //tmp =  xmom_centroid_values[k]*normals[2*ki] + ymom_centroid_values[k]*normals[2*ki+1];
+    //              //if(tmp > 0.){
+    //              //    // Set component of momentum normal to edge to zero
+    //              //    // [equivalently, subtract component of momentum normal to edge]
+    //              //    xmom_centroid_values[k] -= tmp*normals[2*ki]; 
+    //              //    ymom_centroid_values[k] -= tmp*normals[2*ki+1];  
+    //              //    
+    //              //}
+    //           } 
+    //        }
+    //    }
+    //}
+        
+    //   
+    //    
+    //    if((fabs(xmom_explicit_update[k])>1.0e-06)|fabs(ymom_explicit_update[k])>1.0e-06){
+    //        printf(" %e, %e, %e,%e,%e,%e, %e,%e,%e,%d,%e,%e,%e\n", stage_explicit_update[k], xmom_explicit_update[k], ymom_explicit_update[k],
+    //                                    stage_centroid_values[k], bed_centroid_values[k], 
+    //                                    stage_centroid_values[k]- bed_centroid_values[k], 
+    //                                    stage_edge_values[3*neighbours[3*k+0]+neighbour_edges[3*k+0]],
+    //                                    stage_edge_values[3*neighbours[3*k+1]+neighbour_edges[3*k+1]],
+    //                                    stage_edge_values[3*neighbours[3*k+2]+neighbour_edges[3*k+2]],
+    //                                    k,
+    //                                    stage_edge_values[3*k+0],
+    //                                    stage_edge_values[3*k+1],
+    //                                    stage_edge_values[3*k+2]);
+    //    }
+    //}
+    //printf("max ymom edge value is: %e \n", bedtop);
+    //report_python_error(AT, "Written out stage, xmom and ymom updates");
+    //return -1;
+//
+    free(count_wet_edges);
+    //free(count_wet_neighbours);
+    free(edgeflux_store);
+    free(pressuregrad_store);
+
+    return timestep;
+}
+
+// Protect against the water elevation falling below the triangle bed
+int _protect(int N,
+         double minimum_allowed_height,
+         double maximum_allowed_speed,
+         double epsilon,
+         double* wc,
+         double* wv,
+         double* zc,
+         double* zv,
+         double* xmomc,
+         double* ymomc,
+         double* areas,
+         double* xc,
+         double* yc) {
+
+  int k;
+  double hc, bmin, bmax;
+  double u, v, reduced_speed;
+  static double mass_error=0.; 
+  // This acts like minimum_allowed height, but scales with the vertical
+  // distance between the bed_centroid_value and the max bed_edge_value of
+  // every triangle.
+  double minimum_relative_height=0.05; 
+  int mass_added=0;
+
+  // Protect against inifintesimal and negative heights  
+  //if (maximum_allowed_speed < epsilon) {
+    for (k=0; k<N; k++) {
+      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 = 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)) {
+      //  xmomc[k]*=0.1;
+      //  ymomc[k]*=0.1;
+      //}
+
+
+      //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] = 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;
+            //xmomc[k] = xmomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
+            //ymomc[k] = ymomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
+
+        if (hc <= 0.0){
+             // Definine the minimum bed edge value on triangle k.
+             //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]));
+             //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
+             bmin=zc[k]-minimum_allowed_height;//-1.0e-03;
+             // Minimum allowed stage = bmin
+
+             // WARNING: ADDING MASS if wc[k]<bmin
+             if(wc[k]<bmin){
+                 mass_error+=(bmin-wc[k])*areas[k];
+                 mass_added=1; //Flag to warn of added mass                
+
+                 wc[k] = max(wc[k], bmin); 
+                 //printf(" mass_add, %f, %f, %f,%f,%f,%d\n", xc[k], yc[k], wc[k],zc[k],wc[k]-zc[k], k) ;
+             
+
+                 // 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
+                 //wv[3*k] = max(wv[3*k], min(bmin, zv[3*k]-minimum_allowed_height));
+                 //wv[3*k+1] = max(wv[3*k+1], min(bmin, zv[3*k+1]-minimum_allowed_height));
+                 //wv[3*k+2] = max(wv[3*k+2], min(bmin, zv[3*k+2]-minimum_allowed_height));
+                 wv[3*k] = min(bmin, wc[3*k]); //zv[3*k]-minimum_allowed_height);
+                 wv[3*k+1] = min(bmin, wc[3*k+1]); //zv[3*k+1]-minimum_allowed_height);
+                 wv[3*k+2] = min(bmin, wc[3*k+2]); //zv[3*k+2]-minimum_allowed_height);
+            }
+        }
+      }
+    }
+
+  if(mass_added==1){
+     printf("Cumulative mass protection: %f m^3 \n", mass_error);
+  }
+  return 0;
+}
+
+int find_qmin_and_qmax(double dq0, double dq1, double dq2, 
+               double *qmin, double *qmax){
+  // Considering the centroid of an FV triangle and the vertices of its 
+  // auxiliary triangle, find 
+  // qmin=min(q)-qc and qmax=max(q)-qc, 
+  // where min(q) and max(q) are respectively min and max over the
+  // four values (at the centroid of the FV triangle and the auxiliary 
+  // triangle vertices),
+  // and qc is the centroid
+  // dq0=q(vertex0)-q(centroid of FV triangle)
+  // dq1=q(vertex1)-q(vertex0)
+  // dq2=q(vertex2)-q(vertex0)
+
+  // This is a simple implementation 
+  *qmax = max(max(dq0, max(dq0+dq1, dq0+dq2)), 0.0) ;
+  *qmin = min(min(dq0, min(dq0+dq1, dq0+dq2)), 0.0) ;
+ 
+  return 0;
+}
+
+int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
+  // Given provisional jumps dqv from the FV triangle centroid to its 
+  // vertices and jumps qmin (qmax) between the centroid of the FV 
+  // triangle and the minimum (maximum) of the values at the centroid of 
+  // the FV triangle and the auxiliary triangle vertices,
+  // calculate a multiplicative factor phi by which the provisional 
+  // vertex jumps are to be limited
+  
+  int i;
+  double r=1000.0, r0=1.0, phi=1.0;
+  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
+  // FIXME: Perhaps use the epsilon used elsewhere.
+  
+  // Any provisional jump with magnitude < TINY does not contribute to 
+  // the limiting process.
+  //return 0;
+  
+  for (i=0;i<3;i++){
+    if (dqv[i]<-TINY)
+      r0=qmin/dqv[i];
+      
+    if (dqv[i]>TINY)
+      r0=qmax/dqv[i];
+      
+    r=min(r0,r);
+  }
+  
+  phi=min(r*beta_w,1.0);
+  //phi=1.;
+  dqv[0]=dqv[0]*phi;
+  dqv[1]=dqv[1]*phi;
+  dqv[2]=dqv[2]*phi;
+
+  
+  //printf("%lf \n ", beta_w);
+    
+  return 0;
+}
+
+// Computational routine
+int _extrapolate_second_order_edge_sw(int number_of_elements,
+                                 double epsilon,
+                                 double minimum_allowed_height,
+                                 double beta_w,
+                                 double beta_w_dry,
+                                 double beta_uh,
+                                 double beta_uh_dry,
+                                 double beta_vh,
+                                 double beta_vh_dry,
+                                 long* surrogate_neighbours,
+                                 long* neighbour_edges,
+                                 long* number_of_boundaries,
+                                 double* centroid_coordinates,
+                                 double* stage_centroid_values,
+                                 double* xmom_centroid_values,
+                                 double* ymom_centroid_values,
+                                 double* elevation_centroid_values,
+                                 double* height_centroid_values,
+                                 double* edge_coordinates,
+                                 double* stage_edge_values,
+                                 double* xmom_edge_values,
+                                 double* ymom_edge_values,
+                                 double* elevation_edge_values,
+                                 double* height_edge_values,
+                                 double* stage_vertex_values,
+                                 double* xmom_vertex_values,
+                                 double* ymom_vertex_values,
+                                 double* elevation_vertex_values,
+                                 double* height_vertex_values,
+                                 int optimise_dry_cells, 
+                                 int extrapolate_velocity_second_order) {
+                  
+  // Local variables
+  double a, b; // Gradient vector used to calculate edge values from centroids
+  int k, k0, k1, k2, k3, k6, coord_index, i, ii, ktmp, k_wetdry;
+  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
+  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
+  double dqv[3], qmin, qmax, hmin, hmax, bedmax,bedmin, stagemin;
+  double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight, tmp;
+  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vect_norm, l1, l2, a_bs, b_bs, area_e, inv_area_e;
+  
+  double *xmom_centroid_store, *ymom_centroid_store, *min_elevation_edgevalue, *max_elevation_edgevalue;
+  double *cell_wetness_scale;
+  //int *count_wet_neighbours;
+
+  // Use malloc to avoid putting these variables on the stack, which can cause
+  // segfaults in large model runs
+  xmom_centroid_store = malloc(number_of_elements*sizeof(double));
+  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
+  min_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
+  max_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
+  cell_wetness_scale = malloc(number_of_elements*sizeof(double));
+  //count_wet_neighbours = malloc(number_of_elements*sizeof(int));
+ 
+  if(extrapolate_velocity_second_order==1){
+      // Replace momentum centroid with velocity centroid to allow velocity
+      // extrapolation This will be changed back at the end of the routine
+      for (k=0; k<number_of_elements; k++){
+
+          dk = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
+          xmom_centroid_store[k] = xmom_centroid_values[k];
+          xmom_centroid_values[k] = xmom_centroid_values[k]/dk;
+
+          ymom_centroid_store[k] = ymom_centroid_values[k];
+          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]));
+          max_elevation_edgevalue[k] = max(elevation_edge_values[3*k], 
+                                           max(elevation_edge_values[3*k+1],
+                                               elevation_edge_values[3*k+2]));
+
+          // SET HEIGHT IN PLACE
+          height_centroid_values[k] = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
+            
+          }
+
+      }
+
+  //
+  // Compute some useful statistics on wetness/dryness
+  //
+  for(k=0; k<number_of_elements;k++){ 
+      cell_wetness_scale[k] = 0.;
+      // Check if cell k is wet
+      //if(stage_centroid_values[k] > elevation_centroid_values[k]){
+      if(stage_centroid_values[k] > elevation_centroid_values[k] + 2*minimum_allowed_height+epsilon){
+      //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
+          cell_wetness_scale[k] = 1.;  
+      }
+  }
+  //
+  //for(k=0; k<number_of_elements;k++){ 
+  //    // Count the wet neighbours
+  //    count_wet_neighbours[k]=0;
+  //    for (i=0; i<3; i++){
+  //      ktmp = surrogate_neighbours[3*k+i];              
+  //      if( (ktmp >= 0) && cell_wetness_scale[ktmp]>0.){
+  //          count_wet_neighbours[k]+=1;
+  //      }else if(ktmp<=0){
+  //          // Boundary condition -- FIXME: assume wet
+  //          count_wet_neighbours[k]+=1;
+  //      }
+  //       
+  //      
+  //    }
+  //}
+
+  // Alternative 'PROTECT' step
+  for(k=0; k<number_of_elements;k++){ 
+    if(1==1){
+        if((cell_wetness_scale[k]==0. ) ){
+            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
+  //for(k_wetdry=0; k_wetdry<2; k_wetdry++){
+
+      // Begin extrapolation routine
+      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*(1-k_wetdry) ){
+        //  continue;
+        //}
+
+        // Useful indices
+        k3=k*3;
+        k6=k*6;
+
+        if (number_of_boundaries[k]==3)
+        {
+          // No neighbours, set gradient on the triangle to zero
+          
+          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];
+          xmom_edge_values[k3]    = xmom_centroid_values[k];
+          xmom_edge_values[k3+1]  = xmom_centroid_values[k];
+          xmom_edge_values[k3+2]  = xmom_centroid_values[k];
+          ymom_edge_values[k3]    = ymom_centroid_values[k];
+          ymom_edge_values[k3+1]  = ymom_centroid_values[k];
+          ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+          dk = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
+          height_edge_values[k3] = dk;
+          height_edge_values[k3+1] = dk;
+          height_edge_values[k3+2] = dk;
+          
+          continue;
+        }
+        else 
+        {
+          // Triangle k has one or more neighbours. 
+          // Get centroid and edge coordinates of the triangle
+          
+          // Get the edge coordinates
+          xv0 = edge_coordinates[k6];   
+          yv0 = edge_coordinates[k6+1];
+          xv1 = edge_coordinates[k6+2]; 
+          yv1 = edge_coordinates[k6+3];
+          xv2 = edge_coordinates[k6+4]; 
+          yv2 = edge_coordinates[k6+5];
+          
+          // Get the centroid coordinates
+          coord_index = 2*k;
+          x = centroid_coordinates[coord_index];
+          y = centroid_coordinates[coord_index+1];
+          
+          // Store x- and y- differentials for the edges of 
+          // triangle k relative to the centroid
+          dxv0 = xv0 - x; 
+          dxv1 = xv1 - x; 
+          dxv2 = xv2 - x;
+          dyv0 = yv0 - y; 
+          dyv1 = yv1 - y; 
+          dyv2 = yv2 - y;
+
+
+          // Compute bed slope as a reference
+          //area_e = (yv2-yv0)*(xv1-xv0) - (yv1-yv0)*(xv2-xv0);
+          //inv_area_e=1.0/area_e;
+          //a_bs = (yv2 - yv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) -
+          //       (yv1 - yv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
+          //a_bs *= inv_area_e;
+
+          //b_bs = -(xv2 - xv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) +
+          //       (xv1 - xv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
+          //b_bs *= inv_area_e;
+
+          //printf("slopes: %f, %f \n", a_bs, b_bs);
+        }
+
+
+                
+        if (number_of_boundaries[k]<=1)
+        {
+          //==============================================
+          // Number of boundaries <= 1
+          // 'Typical case'
+          //==============================================    
+        
+        
+          // If no boundaries, auxiliary triangle is formed 
+          // from the centroids of the three neighbours
+          // If one boundary, auxiliary triangle is formed 
+          // from this centroid and its two neighbours
+          
+          k0 = surrogate_neighbours[k3];
+          k1 = surrogate_neighbours[k3 + 1];
+          k2 = surrogate_neighbours[k3 + 2];
+         
+          // Test to see whether we accept the surrogate neighbours 
+          // Note that if ki is replaced with k in more than 1 neighbour, then the
+          // triangle area will be zero, and a first order extrapolation will be
+          // 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(k2<0 || 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(k0 < 0 || 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(k1 < 0 || cell_wetness_scale[k1]==0.0){
+              k1 = k ;
+          }
+
+          // Take note if the max neighbour bed elevation is greater than the min
+          // neighbour stage -- suggests a 'steep' bed relative to the flow
+          //bedmax = max(elevation_centroid_values[k], 
+          //             max(elevation_centroid_values[k0],
+          //                 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]))));
+
+          // Get the auxiliary triangle's vertex coordinates 
+          // (really the centroids of neighbouring triangles)
+          coord_index = 2*k0;
+          x0 = centroid_coordinates[coord_index];
+          y0 = centroid_coordinates[coord_index+1];
+          
+          coord_index = 2*k1;
+          x1 = centroid_coordinates[coord_index];
+          y1 = centroid_coordinates[coord_index+1];
+          
+          coord_index = 2*k2;
+          x2 = centroid_coordinates[coord_index];
+          y2 = centroid_coordinates[coord_index+1];
+
+          // Store x- and y- differentials for the vertices 
+          // of the auxiliary triangle
+          dx1 = x1 - x0; 
+          dx2 = x2 - x0;
+          dy1 = y1 - y0; 
+          dy2 = y2 - y0;
+          
+          // Calculate 2*area of the auxiliary triangle
+          // The triangle is guaranteed to be counter-clockwise      
+          area2 = dy2*dx1 - dy1*dx2; 
+           
+          //// Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
+          if ((area2 <= 0.) )// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
+          {
+
+              // Isolated wet cell -- constant stage/depth 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];
+
+              dk=max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
+              height_edge_values[k3] = dk;
+              height_edge_values[k3+1] = dk;
+              height_edge_values[k3+2] = dk;
+
+              xmom_edge_values[k3]    = xmom_centroid_values[k];
+              xmom_edge_values[k3+1]  = xmom_centroid_values[k];
+              xmom_edge_values[k3+2]  = xmom_centroid_values[k];
+              ymom_edge_values[k3]    = ymom_centroid_values[k];
+              ymom_edge_values[k3+1]  = ymom_centroid_values[k];
+              ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+
+              continue;
+          }  
+          
+          // Calculate heights of neighbouring cells
+          hc = stage_centroid_values[k]  - elevation_centroid_values[k];
+          h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
+          h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
+          h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+          hmin = min(min(h0, min(h1, h2)), hc);
+          //// GD FIXME: No longer needed? 
+          //hfactor = 0.0;
+          ////if (hmin > 0.001) 
+          //if (hmin > 0.) 
+          ////if (hc>0.0)
+          //{
+          //  hfactor = 1.0 ;//hmin/(hmin + 0.004);
+            //hfactor=hmin/(hmin + 0.004);
+          //}
+          //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height)), 1.0);
+          hfactor= min(1.5*max(hmin,0.0)/max(hc,minimum_allowed_height), 1.0);
+          
+          //-----------------------------------
+          // stage
+          //-----------------------------------      
+          
+          // Calculate the difference between vertex 0 of the auxiliary 
+          // triangle and the centroid of triangle k
+          dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
+          
+          // Calculate differentials between the vertices 
+          // of the auxiliary triangle (centroids of neighbouring triangles)
+          dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
+          dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
+         
+          inv_area2 = 1.0/area2;
+          // Calculate the gradient of stage on the auxiliary triangle
+          a = dy2*dq1 - dy1*dq2;
+          a *= inv_area2;
+          b = dx1*dq2 - dx2*dq1;
+          b *= inv_area2;
+          // Calculate provisional jumps in stage from the centroid 
+          // of triangle k to its vertices, to be limited
+          dqv[0] = a*dxv0 + b*dyv0;
+          dqv[1] = a*dxv1 + b*dyv1;
+          dqv[2] = a*dxv2 + b*dyv2;
+        
+          // Now we want to find min and max of the centroid and the 
+          // vertices of the auxiliary triangle and compute jumps 
+          // from the centroid to the min and max
+          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+        
+          beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+        
+          
+          // Limit the gradient
+          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];
+
+
+          //-----------------------------------
+          // height
+          //-----------------------------------
+           
+          // Calculate the difference between vertex 0 of the auxiliary 
+          // triangle and the centroid of triangle k
+          dq0 = height_centroid_values[k0] - height_centroid_values[k];
+          
+          // Calculate differentials between the vertices 
+          // of the auxiliary triangle (centroids of neighbouring triangles)
+          dq1 = height_centroid_values[k1] - height_centroid_values[k0];
+          dq2 = height_centroid_values[k2] - height_centroid_values[k0];
+         
+          inv_area2 = 1.0/area2;
+          // Calculate the gradient of height on the auxiliary triangle
+          a = dy2*dq1 - dy1*dq2;
+          a *= inv_area2;
+          b = dx1*dq2 - dx2*dq1;
+          b *= inv_area2;
+          // Calculate provisional jumps in height from the centroid 
+          // of triangle k to its vertices, to be limited
+          dqv[0] = a*dxv0 + b*dyv0;
+          dqv[1] = a*dxv1 + b*dyv1;
+          dqv[2] = a*dxv2 + b*dyv2;
+        
+          // Now we want to find min and max of the centroid and the 
+          // vertices of the auxiliary triangle and compute jumps 
+          // from the centroid to the min and max
+          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+        
+          beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+        
+          
+          // Limit the gradient
+          limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+          height_edge_values[k3+0] = height_centroid_values[k] + dqv[0];
+          height_edge_values[k3+1] = height_centroid_values[k] + dqv[1];
+          height_edge_values[k3+2] = height_centroid_values[k] + dqv[2];
+
+          //-----------------------------------
+          // xmomentum
+          //-----------------------------------            
+
+          // Calculate the difference between vertex 0 of the auxiliary 
+          // triangle and the centroid of triangle k      
+          dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
+          
+          // Calculate differentials between the vertices 
+          // of the auxiliary triangle
+          dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
+          dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
+          
+          // Calculate the gradient of xmom on the auxiliary triangle
+          a = dy2*dq1 - dy1*dq2;
+          a *= inv_area2;
+          b = dx1*dq2 - dx2*dq1;
+          b *= inv_area2;
+          
+          // Calculate provisional jumps in stage from the centroid 
+          // of triangle k to its vertices, to be limited      
+          dqv[0] = a*dxv0+b*dyv0;
+          dqv[1] = a*dxv1+b*dyv1;
+          dqv[2] = a*dxv2+b*dyv2;
+          
+          // Now we want to find min and max of the centroid and the 
+          // vertices of the auxiliary triangle and compute jumps 
+          // from the centroid to the min and max
+          //
+          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+          beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+
+          // Limit the gradient
+          //if((k!=k0)&(k!=k1)&(k!=k2))
+          //if(stagemin>=bedmax){
+          //if((count_wet_neighbours[k]>0)){ // && 
+            // (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{
+          //  dqv[0]=0.;
+          //  dqv[1]=0.;
+          //  dqv[2]=0.;
+          ////  limit_gradient(dqv, qmin, qmax, 0.);
+          //} 
+          //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
+          
+
+          for (i=0; i < 3; i++)
+          {
+            xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
+          }
+          
+          //-----------------------------------
+          // ymomentum
+          //-----------------------------------                  
+
+          // Calculate the difference between vertex 0 of the auxiliary 
+          // triangle and the centroid of triangle k
+          dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
+          
+          // Calculate differentials between the vertices 
+          // of the auxiliary triangle
+          dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
+          dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
+          
+          // Calculate the gradient of xmom on the auxiliary triangle
+          a = dy2*dq1 - dy1*dq2;
+          a *= inv_area2;
+          b = dx1*dq2 - dx2*dq1;
+          b *= inv_area2;
+          
+          // Calculate provisional jumps in stage from the centroid 
+          // of triangle k to its vertices, to be limited
+          dqv[0] = a*dxv0 + b*dyv0;
+          dqv[1] = a*dxv1 + b*dyv1;
+          dqv[2] = a*dxv2 + b*dyv2;
+          
+          // Now we want to find min and max of the centroid and the 
+          // vertices of the auxiliary triangle and compute jumps 
+          // from the centroid to the min and max
+          //
+          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+          
+          beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
+
+          // Limit the gradient
+          //if(stagemin>=bedmax){
+          //if((count_wet_neighbours[k]>0)){//&&
+             //(cell_wetness_scale[k]==1.0)){
+            //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
+          //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
+            limit_gradient(dqv, qmin, qmax, beta_tmp);
+          //}else{
+          //  dqv[0]=0.;
+          //  dqv[1]=0.;
+          //  dqv[2]=0.;
+          //} 
+          
+          for (i=0;i<3;i++)
+          {
+            ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+          }
+          
+        } // End number_of_boundaries <=1 
+        else
+        {
+
+          //==============================================
+          // Number of boundaries == 2
+          //==============================================        
+            
+          // One internal neighbour and gradient is in direction of the neighbour's centroid
+          
+          // Find the only internal neighbour (k1?)
+          for (k2 = k3; k2 < k3 + 3; k2++)
+          {
+          // Find internal neighbour of triangle k      
+          // k2 indexes the edges of triangle k   
+          
+              if (surrogate_neighbours[k2] != k)
+              {
+                 break;
+              }
+          }
+          
+          if ((k2 == k3 + 3)) 
+          {
+            // If we didn't find an internal neighbour
+            report_python_error(AT, "Internal neighbour not found");
+            return -1;
+          }
+          
+          k1 = surrogate_neighbours[k2];
+          
+          // The coordinates of the triangle are already (x,y). 
+          // Get centroid of the neighbour (x1,y1)
+          coord_index = 2*k1;
+          x1 = centroid_coordinates[coord_index];
+          y1 = centroid_coordinates[coord_index + 1];
+          
+          // Compute x- and y- distances between the centroid of 
+          // triangle k and that of its neighbour
+          dx1 = x1 - x; 
+          dy1 = y1 - y;
+          
+          // Set area2 as the square of the distance
+          area2 = dx1*dx1 + dy1*dy1;
+          
+          // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) 
+          // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
+          // respectively correspond to the x- and y- gradients 
+          // of the conserved quantities
+          dx2 = 1.0/area2;
+          dy2 = dx2*dy1;
+          dx2 *= dx1;
+          
+          
+          //-----------------------------------
+          // stage
+          //-----------------------------------            
+
+          // Compute differentials
+          dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
+          
+          // Calculate the gradient between the centroid of triangle k 
+          // and that of its neighbour
+          a = dq1*dx2;
+          b = dq1*dy2;
+          
+          // Calculate provisional edge jumps, to be limited
+          dqv[0] = a*dxv0 + b*dyv0;
+          dqv[1] = a*dxv1 + b*dyv1;
+          dqv[2] = a*dxv2 + b*dyv2;
+          
+          // Now limit the jumps
+          if (dq1>=0.0)
+          {
+            qmin=0.0;
+            qmax=dq1;
+          }
+          else
+          {
+            qmin = dq1;
+            qmax = 0.0;
+          }
+          
+          // Limit the gradient
+          limit_gradient(dqv, qmin, qmax, beta_w);
+          
+          stage_edge_values[k3] = 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];
+
+          //-----------------------------------
+          // height
+          //-----------------------------------
+          
+          // Compute differentials
+          dq1 = height_centroid_values[k1] - height_centroid_values[k];
+          
+          // Calculate the gradient between the centroid of triangle k 
+          // and that of its neighbour
+          a = dq1*dx2;
+          b = dq1*dy2;
+          
+          // Calculate provisional edge jumps, to be limited
+          dqv[0] = a*dxv0 + b*dyv0;
+          dqv[1] = a*dxv1 + b*dyv1;
+          dqv[2] = a*dxv2 + b*dyv2;
+          
+          // Now limit the jumps
+          if (dq1>=0.0)
+          {
+            qmin=0.0;
+            qmax=dq1;
+          }
+          else
+          {
+            qmin = dq1;
+            qmax = 0.0;
+          }
+          
+          // Limit the gradient
+          limit_gradient(dqv, qmin, qmax, beta_w);
+          
+          height_edge_values[k3] = height_centroid_values[k] + dqv[0];
+          height_edge_values[k3 + 1] = height_centroid_values[k] + dqv[1];
+          height_edge_values[k3 + 2] = height_centroid_values[k] + dqv[2];
+
+          //-----------------------------------
+          // xmomentum
+          //-----------------------------------                        
+          
+          // Compute differentials
+          dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
+          
+          // Calculate the gradient between the centroid of triangle k 
+          // and that of its neighbour
+          a = dq1*dx2;
+          b = dq1*dy2;
+          
+          // Calculate provisional edge jumps, to be limited
+          dqv[0] = a*dxv0+b*dyv0;
+          dqv[1] = a*dxv1+b*dyv1;
+          dqv[2] = a*dxv2+b*dyv2;
+          
+          // Now limit the jumps
+          if (dq1 >= 0.0)
+          {
+            qmin = 0.0;
+            qmax = dq1;
+          }
+          else
+          {
+            qmin = dq1;
+            qmax = 0.0;
+          }
+          
+          // Limit the gradient      
+          limit_gradient(dqv, qmin, qmax, beta_w);
+          
+          //for (i=0;i<3;i++)
+          //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
+          //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
+          //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
+          
+          for (i = 0; i < 3;i++)
+          {
+              xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
+          }
+          
+          //-----------------------------------
+          // ymomentum
+          //-----------------------------------                        
+
+          // Compute differentials
+          dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
+          
+          // Calculate the gradient between the centroid of triangle k 
+          // and that of its neighbour
+          a = dq1*dx2;
+          b = dq1*dy2;
+          
+          // Calculate provisional edge jumps, to be limited
+          dqv[0] = a*dxv0 + b*dyv0;
+          dqv[1] = a*dxv1 + b*dyv1;
+          dqv[2] = a*dxv2 + b*dyv2;
+          
+          // Now limit the jumps
+          if (dq1>=0.0)
+          {
+            qmin = 0.0;
+            qmax = dq1;
+          } 
+          else
+          {
+            qmin = dq1;
+            qmax = 0.0;
+          }
+          
+          // Limit the gradient
+          limit_gradient(dqv, qmin, qmax, beta_w);
+          
+          for (i=0;i<3;i++)
+                  {
+                  ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+                  }
+        } // else [number_of_boundaries==2]
+      } // for k=0 to number_of_elements-1
+  //} 
+
+  //for(k=0;k<number_of_elements; k++){
+  //    // 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(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(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]);
+  //          }
+  //      }
+  //    }
+  //}
+
+  // Compute vertex values of quantities
+  for (k=0; k<number_of_elements; k++){
+      k3=3*k;
+
+      // Compute stage vertex values 
+      stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ; 
+      stage_vertex_values[k3+1] =  stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1]; 
+      stage_vertex_values[k3+2] =  stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2]; 
+      
+      // Compute height vertex values 
+      height_vertex_values[k3] = height_edge_values[k3+1] + height_edge_values[k3+2] -height_edge_values[k3] ; 
+      height_vertex_values[k3+1] =  height_edge_values[k3] + height_edge_values[k3+2]-height_edge_values[k3+1]; 
+      height_vertex_values[k3+2] =  height_edge_values[k3] + height_edge_values[k3+1]-height_edge_values[k3+2]; 
+
+     // Compute xmom vertex values 
+      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ; 
+      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1]; 
+      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2]; 
+
+      // Compute ymom vertex values 
+      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ; 
+      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1]; 
+      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2]; 
+
+      // If needed, convert from velocity to momenta
+      if(extrapolate_velocity_second_order==1){
+          //Convert velocity back to momenta at centroids
+          xmom_centroid_values[k] = xmom_centroid_store[k];
+          ymom_centroid_values[k] = ymom_centroid_store[k];
+
+          // Re-compute momenta at edges
+          for (i=0; i<3; i++){
+              de[i] = height_edge_values[k3+i]; //max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
+              //if(de[i]<minimum_allowed_height){
+              //  de[i]=0.;
+              //}
+              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
+              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i];
+          }
+
+          // Re-compute momenta at vertices
+          for (i=0; i<3; i++){
+              de[i] = height_vertex_values[k3+i]; //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
+              xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
+              ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
+          }
+
+      }
+
+      // Compute new bed elevation
+      elevation_edge_values[k3]=stage_edge_values[k3]-height_edge_values[k3];
+      elevation_edge_values[k3+1]=stage_edge_values[k3+1]-height_edge_values[k3+1];
+      elevation_edge_values[k3+2]=stage_edge_values[k3+2]-height_edge_values[k3+2];
+      elevation_vertex_values[k3]=stage_vertex_values[k3]-height_vertex_values[k3];
+      elevation_vertex_values[k3+1]=stage_vertex_values[k3+1]-height_vertex_values[k3+1];
+      elevation_vertex_values[k3+2]=stage_vertex_values[k3+2]-height_vertex_values[k3+2];
+
+   
+  } 
+
+  free(xmom_centroid_store);
+  free(ymom_centroid_store);
+  free(min_elevation_edgevalue);
+  free(max_elevation_edgevalue);
+  free(cell_wetness_scale);
+  //free(count_wet_neighbours);
+
+  return 0;
+}           
+
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+
+//========================================================================
+// Compute fluxes
+//========================================================================
+
+// Modified central flux function
+
+PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function_central"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Python call:
+    domain.timestep = compute_fluxes(timestep,
+                                     domain.epsilon,
+                     domain.H0,
+                                     domain.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,
+                                     domain.radii,
+                                     domain.areas,
+                                     tri_full_flag,
+                                     Stage.edge_values,
+                                     Xmom.edge_values,
+                                     Ymom.edge_values,
+                                     Bed.edge_values,
+                                     Stage.boundary_values,
+                                     Xmom.boundary_values,
+                                     Ymom.boundary_values,
+                                     Stage.explicit_update,
+                                     Xmom.explicit_update,
+                                     Ymom.explicit_update,
+                                     already_computed_flux,
+                                     optimise_dry_cells,
+                                     stage.centroid_values,
+                                     bed.centroid_values)
+
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+
+
+  */
+
+
+  PyArrayObject *boundary_flux_sum, *neighbours, *neighbour_edges,
+    *normals, *edgelengths, *radii, *areas,
+    *tri_full_flag,
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *bed_edge_values,
+    *height_edge_values,
+    *stage_boundary_values,
+    *xmom_boundary_values,
+    *ymom_boundary_values,
+    *boundary_flux_type,
+    *stage_explicit_update,
+    *xmom_explicit_update,
+    *ymom_explicit_update,
+    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
+    *max_speed_array, //Keeps track of max speeds for each triangle
+    *stage_centroid_values,
+    *xmom_centroid_values,
+    *ymom_centroid_values,
+    *bed_centroid_values,
+    *height_centroid_values,
+    *bed_vertex_values;
+    
+  double timestep, epsilon, H0, g;
+  int optimise_dry_cells, timestep_order;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOOiiOOOOOO",
+            &timestep,
+            &epsilon,
+            &H0,
+            &g,
+            &boundary_flux_sum,
+            &neighbours,
+            &neighbour_edges,
+            &normals,
+            &edgelengths, &radii, &areas,
+            &tri_full_flag,
+            &stage_edge_values,
+            &xmom_edge_values,
+            &ymom_edge_values,
+            &bed_edge_values,
+            &height_edge_values,
+            &stage_boundary_values,
+            &xmom_boundary_values,
+            &ymom_boundary_values,
+            &boundary_flux_type,
+            &stage_explicit_update,
+            &xmom_explicit_update,
+            &ymom_explicit_update,
+            &already_computed_flux,
+            &max_speed_array,
+            &optimise_dry_cells,
+            &timestep_order,
+            &stage_centroid_values,
+            &xmom_centroid_values,
+            &ymom_centroid_values,
+            &bed_centroid_values,
+            &height_centroid_values,
+            &bed_vertex_values)) {
+    report_python_error(AT, "could not parse input arguments");
+    return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(neighbours);
+  CHECK_C_CONTIG(neighbour_edges);
+  CHECK_C_CONTIG(normals);
+  CHECK_C_CONTIG(edgelengths);
+  CHECK_C_CONTIG(radii);
+  CHECK_C_CONTIG(areas);
+  CHECK_C_CONTIG(tri_full_flag);
+  CHECK_C_CONTIG(stage_edge_values);
+  CHECK_C_CONTIG(xmom_edge_values);
+  CHECK_C_CONTIG(ymom_edge_values);
+  CHECK_C_CONTIG(bed_edge_values);
+  CHECK_C_CONTIG(height_edge_values);
+  CHECK_C_CONTIG(stage_boundary_values);
+  CHECK_C_CONTIG(xmom_boundary_values);
+  CHECK_C_CONTIG(ymom_boundary_values);
+  CHECK_C_CONTIG(boundary_flux_type);
+  CHECK_C_CONTIG(stage_explicit_update);
+  CHECK_C_CONTIG(xmom_explicit_update);
+  CHECK_C_CONTIG(ymom_explicit_update);
+  CHECK_C_CONTIG(already_computed_flux);
+  CHECK_C_CONTIG(max_speed_array);
+  CHECK_C_CONTIG(stage_centroid_values);
+  CHECK_C_CONTIG(xmom_centroid_values);
+  CHECK_C_CONTIG(ymom_centroid_values);
+  CHECK_C_CONTIG(bed_centroid_values);
+  CHECK_C_CONTIG(height_centroid_values);
+  CHECK_C_CONTIG(bed_vertex_values);
+
+  int number_of_elements = stage_edge_values -> dimensions[0];
+
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  timestep = _compute_fluxes_central(number_of_elements,
+                     timestep,
+                     epsilon,
+                     H0,
+                     g,
+                     (double*) boundary_flux_sum -> data,
+                     (long*) neighbours -> data,
+                     (long*) neighbour_edges -> data,
+                     (double*) normals -> data,
+                     (double*) edgelengths -> data, 
+                     (double*) radii -> data, 
+                     (double*) areas -> data,
+                     (long*) tri_full_flag -> data,
+                     (double*) stage_edge_values -> data,
+                     (double*) xmom_edge_values -> data,
+                     (double*) ymom_edge_values -> data,
+                     (double*) bed_edge_values -> data,
+                     (double*) height_edge_values -> data,
+                     (double*) stage_boundary_values -> data,
+                     (double*) xmom_boundary_values -> data,
+                     (double*) ymom_boundary_values -> data,
+                     (long*)   boundary_flux_type -> data,
+                     (double*) stage_explicit_update -> data,
+                     (double*) xmom_explicit_update -> data,
+                     (double*) ymom_explicit_update -> data,
+                     (long*) already_computed_flux -> data,
+                     (double*) max_speed_array -> data,
+                     optimise_dry_cells, 
+                     timestep_order,
+                     (double*) stage_centroid_values -> data, 
+                     (double*) xmom_centroid_values -> data, 
+                     (double*) ymom_centroid_values -> data, 
+                     (double*) bed_centroid_values -> data,
+                     (double*) height_centroid_values -> data,
+                     (double*) bed_vertex_values -> data);
+  
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
+PyObject *flux_function_central(PyObject *self, PyObject *args) {
+  //
+  // Gateway to innermost flux function.
+  // This is only used by the unit tests as the c implementation is
+  // normally called by compute_fluxes in this module.
+
+
+  PyArrayObject *normal, *ql, *qr,  *edgeflux;
+  double g, epsilon, max_speed, H0, zl, zr;
+  double h0, limiting_threshold, pressure_flux, smooth;
+  double hle, hre;
+
+  if (!PyArg_ParseTuple(args, "OOOddOddd",
+            &normal, &ql, &qr, &zl, &zr, &edgeflux,
+            &epsilon, &g, &H0)) {
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
+  }
+
+  printf("Warning: This interface is broken -- do not use \n");
+  
+  h0 = H0; 
+  hle=1.0; // Fake values to force this to compile
+  hre=1.0; // Fake values to force this to compile	      
+  limiting_threshold = 10*H0; // Avoid applying limiter below this
+                              // threshold for performance reasons.
+                              // See ANUGA manual under flux limiting  
+ 
+  pressure_flux = 0.0; // Store the water pressure related component of the flux 
+  smooth = 1.0 ; // term to scale diffusion in riemann solver
+
+  _flux_function_central((double*) ql -> data, 
+			 (double*) qr -> data, 
+			 zl, 
+			 zr,
+             hle,
+             hre,                         
+			 ((double*) normal -> data)[0],
+			 ((double*) normal -> data)[1],          
+			 epsilon, h0, limiting_threshold,
+			 g,
+			 (double*) edgeflux -> data, 
+			 &max_speed,
+             &pressure_flux,
+			 ((double*) normal -> data)[0],
+			 ((double*) normal -> data)[1], 
+			 ((double*) normal -> data)[1]
+             );
+  
+  return Py_BuildValue("d", max_speed);  
+}
+
+//========================================================================
+// Gravity
+//========================================================================
+
+PyObject *gravity(PyObject *self, PyObject *args) {
+  //
+  //  gravity(g, h, v, x, xmom, ymom)
+  //
+  
+  
+  PyArrayObject *h, *v, *x, *xmom, *ymom;
+  int k, N, k3, k6;
+  double g, avg_h, zx, zy;
+  double x0, y0, x1, y1, x2, y2, z0, z1, z2;
+  //double epsilon;
+  
+  if (!PyArg_ParseTuple(args, "dOOOOO",
+			&g, &h, &v, &x,
+			&xmom, &ymom)) {
+    //&epsilon)) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
+    return NULL;
+  }
+  
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(h);
+  CHECK_C_CONTIG(v);
+  CHECK_C_CONTIG(x);
+  CHECK_C_CONTIG(xmom);
+  CHECK_C_CONTIG(ymom);
+  
+  N = h -> dimensions[0];
+  for (k=0; k<N; k++) {
+    k3 = 3*k;  // base index
+    
+    // Get bathymetry
+    z0 = ((double*) v -> data)[k3 + 0];
+    z1 = ((double*) v -> data)[k3 + 1];
+    z2 = ((double*) v -> data)[k3 + 2];
+    
+    // Optimise for flat bed
+    // Note (Ole): This didn't produce measurable speed up.
+    // Revisit later
+    //if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
+    //  continue;
+    //} 
+    
+    // Get average depth from centroid values
+    avg_h = ((double *) h -> data)[k];
+    
+    // Compute bed slope
+    k6 = 6*k;  // base index
+    
+    x0 = ((double*) x -> data)[k6 + 0];
+    y0 = ((double*) x -> data)[k6 + 1];
+    x1 = ((double*) x -> data)[k6 + 2];
+    y1 = ((double*) x -> data)[k6 + 3];
+    x2 = ((double*) x -> data)[k6 + 4];
+    y2 = ((double*) x -> data)[k6 + 5];
+    
+    
+    _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
+    
+    // Update momentum
+    ((double*) xmom -> data)[k] += -g*zx*avg_h;
+    ((double*) ymom -> data)[k] += -g*zy*avg_h;
+  }
+  
+  return Py_BuildValue("");
+}
+
+
+PyObject *extrapolate_second_order_edge_sw(PyObject *self, PyObject *args) {
+  /*Compute the edge values based on a linear reconstruction 
+    on each triangle
+    
+    Python call:
+    extrapolate_second_order_sw(domain.surrogate_neighbours,
+                                domain.number_of_boundaries
+                                domain.centroid_coordinates,
+                                Stage.centroid_values
+                                Xmom.centroid_values
+                                Ymom.centroid_values
+                                domain.edge_coordinates,
+                                Stage.edge_values,
+                                Xmom.edge_values,
+                                Ymom.edge_values)
+
+    Post conditions:
+        The edges of each triangle have values from a 
+        limited linear reconstruction
+        based on centroid values
+
+  */
+  PyArrayObject *surrogate_neighbours,
+    *neighbour_edges,
+    *number_of_boundaries,
+    *centroid_coordinates,
+    *stage_centroid_values,
+    *xmom_centroid_values,
+    *ymom_centroid_values,
+    *elevation_centroid_values,
+    *height_centroid_values,
+    *edge_coordinates,
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *elevation_edge_values,
+    *height_edge_values,
+    *stage_vertex_values,
+    *xmom_vertex_values,
+    *ymom_vertex_values,
+    *elevation_vertex_values,
+    *height_vertex_values;
+  
+  PyObject *domain;
+
+  
+  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;    
+  double minimum_allowed_height, epsilon;
+  int optimise_dry_cells, number_of_elements, extrapolate_velocity_second_order, e, e2;
+  
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOOOOii",
+            &domain,
+            &surrogate_neighbours,
+            &neighbour_edges,
+            &number_of_boundaries,
+            &centroid_coordinates,
+            &stage_centroid_values,
+            &xmom_centroid_values,
+            &ymom_centroid_values,
+            &elevation_centroid_values,
+            &height_centroid_values,
+            &edge_coordinates,
+            &stage_edge_values,
+            &xmom_edge_values,
+            &ymom_edge_values,
+            &elevation_edge_values,
+            &height_edge_values,
+            &stage_vertex_values,
+            &xmom_vertex_values,
+            &ymom_vertex_values,
+            &elevation_vertex_values,
+            &height_vertex_values,
+            &optimise_dry_cells,
+            &extrapolate_velocity_second_order)) {         
+
+    report_python_error(AT, "could not parse input arguments");
+    return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(surrogate_neighbours);
+  CHECK_C_CONTIG(neighbour_edges);
+  CHECK_C_CONTIG(number_of_boundaries);
+  CHECK_C_CONTIG(centroid_coordinates);
+  CHECK_C_CONTIG(stage_centroid_values);
+  CHECK_C_CONTIG(xmom_centroid_values);
+  CHECK_C_CONTIG(ymom_centroid_values);
+  CHECK_C_CONTIG(elevation_centroid_values);
+  CHECK_C_CONTIG(height_centroid_values);
+  CHECK_C_CONTIG(edge_coordinates);
+  CHECK_C_CONTIG(stage_edge_values);
+  CHECK_C_CONTIG(xmom_edge_values);
+  CHECK_C_CONTIG(ymom_edge_values);
+  CHECK_C_CONTIG(elevation_edge_values);
+  CHECK_C_CONTIG(height_edge_values);
+  CHECK_C_CONTIG(stage_vertex_values);
+  CHECK_C_CONTIG(xmom_vertex_values);
+  CHECK_C_CONTIG(ymom_vertex_values);
+  CHECK_C_CONTIG(elevation_vertex_values);
+  CHECK_C_CONTIG(height_vertex_values);
+  
+  // Get the safety factor beta_w, set in the config.py file. 
+  // This is used in the limiting process
+  
+
+  beta_w                 = get_python_double(domain,"beta_w");
+  beta_w_dry             = get_python_double(domain,"beta_w_dry");
+  beta_uh                = get_python_double(domain,"beta_uh");
+  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
+  beta_vh                = get_python_double(domain,"beta_vh");
+  beta_vh_dry            = get_python_double(domain,"beta_vh_dry");  
+
+  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
+  epsilon                = get_python_double(domain,"epsilon");
+
+  number_of_elements = stage_centroid_values -> dimensions[0];  
+
+  //printf("In C before Extrapolate");
+  //e=1;
+  // Call underlying computational routine
+  e = _extrapolate_second_order_edge_sw(number_of_elements,
+                   epsilon,
+                   minimum_allowed_height,
+                   beta_w,
+                   beta_w_dry,
+                   beta_uh,
+                   beta_uh_dry,
+                   beta_vh,
+                   beta_vh_dry,
+                   (long*) surrogate_neighbours -> data,
+                   (long*) neighbour_edges -> data,
+                   (long*) number_of_boundaries -> data,
+                   (double*) centroid_coordinates -> data,
+                   (double*) stage_centroid_values -> data,
+                   (double*) xmom_centroid_values -> data,
+                   (double*) ymom_centroid_values -> data,
+                   (double*) elevation_centroid_values -> data,
+                   (double*) height_centroid_values -> data,
+                   (double*) edge_coordinates -> data,
+                   (double*) stage_edge_values -> data,
+                   (double*) xmom_edge_values -> data,
+                   (double*) ymom_edge_values -> data,
+                   (double*) elevation_edge_values -> data,
+                   (double*) height_edge_values -> data,
+                   (double*) stage_vertex_values -> data,
+                   (double*) xmom_vertex_values -> data,
+                   (double*) ymom_vertex_values -> data,
+                   (double*) elevation_vertex_values -> data,
+                   (double*) height_vertex_values -> data,
+                   optimise_dry_cells, 
+                   extrapolate_velocity_second_order);
+
+  if (e == -1) {
+    // Use error string set inside computational routine
+    return NULL;
+  }                   
+  
+  
+  return Py_BuildValue("");
+  
+}// extrapolate_second-order_edge_sw
+//========================================================================
+// Protect -- to prevent the water level from falling below the minimum 
+// bed_edge_value
+//========================================================================
+
+PyObject *protect(PyObject *self, PyObject *args) {
+  //
+  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
+
+
+  PyArrayObject
+  *wc,            //Stage at centroids
+  *wv,            //Stage at vertices
+  *zc,            //Elevation at centroids
+  *zv,            //Elevation at vertices
+  *xmomc,         //Momentums at centroids
+  *ymomc,
+  *areas,         // Triangle areas
+  *xc,
+  *yc;
+
+  int N;
+  double minimum_allowed_height, maximum_allowed_speed, epsilon;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddOOOOOOOOO",
+            &minimum_allowed_height,
+            &maximum_allowed_speed,
+            &epsilon,
+            &wc, &wv, &zc, &zv, &xmomc, &ymomc, &areas, &xc, &yc)) {
+    report_python_error(AT, "could not parse input arguments");
+    return NULL;
+  }  
+
+  N = wc -> dimensions[0];
+
+  _protect(N,
+       minimum_allowed_height,
+       maximum_allowed_speed,
+       epsilon,
+       (double*) wc -> data,
+       (double*) wv -> data,
+       (double*) zc -> data,
+       (double*) zv -> data,
+       (double*) xmomc -> data,
+       (double*) ymomc -> data,
+       (double*) areas -> data,
+       (double*) xc -> data,
+       (double*) yc -> data );
+
+  return Py_BuildValue("");
+}
+
+//========================================================================
+// Method table for python module
+//========================================================================
+
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
+  {"gravity_c",        gravity,            METH_VARARGS, "Print out"},
+  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},  
+  {"extrapolate_second_order_edge_sw", extrapolate_second_order_edge_sw, METH_VARARGS, "Print out"},
+  {"protect",          protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {NULL, NULL, 0, NULL}
+};
+
+// Module initialisation
+void initswb2_domain_ext(void){
+  Py_InitModule("swb2_domain_ext", MethodTable);
+
+  import_array(); // Necessary for handling of NumPY structures
+}
Index: trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py	(revision 8990)
@@ -14,5 +14,6 @@
 #from anuga import *
 #from balanced_basic import *
-from balanced_dev import *
+#from balanced_dev import *
+from bal_and import *
 #from anuga_tsunami import *
 #from balanced_basic.swb2_domain import *
Index: trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 8990)
@@ -1,8 +1,10 @@
-import util
+#import util
+from bal_and import plot_utils as util
 from matplotlib import pyplot as pyplot
 #import pylab
 
 # Time-index to plot outputs from
-index=900
+index=150
+#index=900
 
 #p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
@@ -28,5 +30,5 @@
 
 # Numerical results along a central channel 'slice'
-V1 = p.stage[index,v] - p.elev[v]
+V1 = p.stage[index,v] - p.elev[index,v]
 V2 = p.yvel[index,v]
 V3 = p.xvel[index,v]
@@ -90,11 +92,11 @@
 
 pyplot.clf()
-analytical_stage = min(p.elev[v1]) + dc_analytical
-analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
-analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+analytical_stage = min(p.elev[index,v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[index,v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[index,v1])
 pyplot.figure(figsize=(12.,8.))
 pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
 pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.elev[index,v1],'o', label='bed elevation (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
@@ -109,11 +111,11 @@
 
 pyplot.clf()
-analytical_stage = min(p.elev[v1]) + dc_analytical
-analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
-analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+analytical_stage = min(p.elev[index,v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[index,v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[index,v1])
 pyplot.figure(figsize=(12.,8.))
 pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
 pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.elev[index,v1],'o', label='bed elevation (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
@@ -127,11 +129,11 @@
 
 pyplot.clf()
-analytical_stage = min(p.elev[v1]) + dc_analytical
-analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
-analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+analytical_stage = min(p.elev[index,v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[index,v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[index,v1])
 pyplot.figure(figsize=(12.,8.))
 pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
 pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.elev[index,v1],'o', label='bed elevation (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
Index: trunk/anuga_work/development/gareth/tests/parabolic/parabolic.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/parabolic/parabolic.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/parabolic/parabolic.py	(revision 8990)
@@ -8,9 +8,11 @@
 import numpy
 
+from bal_and import *
+
 #from balanced_dev import *
 #from anuga_tsunami import *
 #from balanced_basic import *
 #from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
-from anuga.shallow_water.shallow_water_domain import Domain as Domain
+#from anuga.shallow_water.shallow_water_domain import Domain as Domain
 #---------
 #Setup computational domain
@@ -21,7 +23,7 @@
 domain.set_name('parabola_v2')                         # Output to file runup.sww
 domain.set_datadir('.')                          # Use current folder
-domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
 domain.set_minimum_allowed_height(0.001)
-domain.set_flow_algorithm('tsunami')
+
+#domain.set_flow_algorithm('tsunami')
 # Time stepping
 #domain.set_timestepping_method('euler') # Default
Index: trunk/anuga_work/development/gareth/tests/runup/runup.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/runup/runup.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/runup/runup.py	(revision 8990)
@@ -18,5 +18,5 @@
 #path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
 #from balanced_basic import *
-from balanced_dev import *
+from bal_and import *
 #from anuga_tsunami import *
 #---------
@@ -28,5 +28,5 @@
 domain.set_name('runup_v2')                         # Output to file runup.sww
 domain.set_datadir('.')                          # Use current folder
-domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
+#domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
 #domain.set_store_vertices_uniquely(True)
 #------------------
@@ -34,5 +34,5 @@
 #------------------
 
-scale_me=100.0
+scale_me=1.0
 
 def topography(x,y):
@@ -57,5 +57,5 @@
 Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
 Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
-Bd=anuga.Dirichlet_boundary([-0.2*scale_me,0.,0.])       # Constant boundary values -- not used in this example
+Bd=anuga.Dirichlet_boundary([-0.18*scale_me,0.,0.])       # Constant boundary values -- not used in this example
 #Bw=anuga.Time_boundary(domain=domain,
 #	f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
@@ -64,5 +64,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})
 
 #------------------------------
@@ -70,9 +70,9 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.2,finaltime=60.0):
+for t in domain.evolve(yieldstep=0.2,finaltime=10.0):
     print domain.timestepping_statistics()
     uh=domain.quantities['xmomentum'].centroid_values
     vh=domain.quantities['ymomentum'].centroid_values
-    depth=domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
+    depth=domain.quantities['height'].centroid_values#domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
     depth=depth*(depth>1.0e-06) + 1.0e-06
     vel=((uh/depth)**2 + (vh/depth)**2)**0.5
Index: trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 8990)
@@ -14,5 +14,6 @@
 #from swb2_domain import *
 #from balanced_basic import *
-from balanced_dev import *
+#from balanced_dev import *
+from bal_and import *
 #from anuga_tsunami import *
 
Index: trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py	(revision 8990)
@@ -4,5 +4,6 @@
 # Import Modules
 #---------------
-import util
+#import util
+from bal_and 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_sparse.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py	(revision 8990)
@@ -21,5 +21,5 @@
 from anuga.structures.inlet_operator import Inlet_operator
 from anuga.shallow_water.shallow_water_domain import Domain as Domain
-from balanced_dev import *
+from bal_and import *
 #from balanced_dev import Domain as Domain
 #from anuga_tsunami import *
Index: trunk/anuga_work/development/gareth/tests/wave/run_wave.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/wave/run_wave.py	(revision 8893)
+++ trunk/anuga_work/development/gareth/tests/wave/run_wave.py	(revision 8990)
@@ -15,5 +15,6 @@
 
 #from balanced_dev import *
-from anuga_tsunami import *
+#from anuga_tsunami import *
+from bal_and import *
 #from balanced_dev import Domain as Domain
 
@@ -89,7 +90,7 @@
 amplitude = 1
 wave_length = 300.0
-Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary
-## Sine wave
-                  f=lambda t: [(-amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
+#Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary
+### Sine wave
+#                  f=lambda t: [(-amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
 ## Sawtooth?
 #                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
@@ -129,5 +130,5 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep = 10.0, finaltime = 60*60.*3.):
+for t in domain.evolve(yieldstep = 10.0, finaltime = 5*10*60.*3.):
     domain.write_time()
     if interactive_visualisation:
