Index: /trunk/anuga_work/development/gareth/balanced_basic/__init__.py
===================================================================
--- /trunk/anuga_work/development/gareth/balanced_basic/__init__.py	(revision 8302)
+++ /trunk/anuga_work/development/gareth/balanced_basic/__init__.py	(revision 8302)
@@ -0,0 +1,23 @@
+"""Make directory available as a Python package
+"""
+
+# Add path of package to PYTHONPATH to allow C-extensions to be loaded
+import sys
+sys.path += __path__
+
+# Make selected classes available directly
+from swb2_domain import Domain
+#, \
+#     Transmissive_boundary,\
+#     Reflective_boundary,\
+#     Dirichlet_boundary, Time_boundary, File_boundary,\
+#     Transmissive_momentum_set_stage_boundary,\
+#     Dirichlet_discharge_boundary,\
+#     Field_boundary,\
+#     Transmissive_stage_zero_momentum_boundary,\
+#     Transmissive_n_momentum_zero_t_momentum_set_stage_boundary,\
+#     AWI_boundary
+#
+
+
+ 
Index: /trunk/anuga_work/development/gareth/balanced_basic/compile_all.py
===================================================================
--- /trunk/anuga_work/development/gareth/balanced_basic/compile_all.py	(revision 8302)
+++ /trunk/anuga_work/development/gareth/balanced_basic/compile_all.py	(revision 8302)
@@ -0,0 +1,23 @@
+import os
+import time
+
+buildroot = '..' + os.sep + '..' + os.sep + '..' + os.sep + '..' + os.sep + 'anuga_core'
+
+compile_script = buildroot + os.sep + 'source' + os.sep + 'anuga' + \
+                 os.sep + 'utilities' + os.sep + 'compile.py'
+
+t0 = time.time()
+
+# Attempt to compile all ANUGA extensions
+
+execfile(compile_script)
+
+
+  
+
+
+print        
+print 'That took %.3fs' %(time.time() - t0)
+
+if sys.platform == 'win32':
+    raw_input('Press the RETURN key')
Index: /trunk/anuga_work/development/gareth/balanced_basic/swb2_boundary_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/balanced_basic/swb2_boundary_conditions.py	(revision 8302)
+++ /trunk/anuga_work/development/gareth/balanced_basic/swb2_boundary_conditions.py	(revision 8302)
@@ -0,0 +1,33 @@
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Boundary
+
+class Transmissive_boundary(Boundary):
+    """Transmissive boundary returns same conserved quantities as
+    those present in its neighbour volume.
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain = None):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for transmissive boundary'
+            raise Exception, msg
+
+        self.domain = domain
+
+    def __repr__(self):
+        return 'Transmissive_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
Index: /trunk/anuga_work/development/gareth/balanced_basic/swb2_domain.py
===================================================================
--- /trunk/anuga_work/development/gareth/balanced_basic/swb2_domain.py	(revision 8302)
+++ /trunk/anuga_work/development/gareth/balanced_basic/swb2_domain.py	(revision 8302)
@@ -0,0 +1,390 @@
+
+from anuga.shallow_water.shallow_water_domain import *
+from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain
+
+
+##############################################################################
+# 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)
+        
+        #---------------------
+        # set some defaults
+        #---------------------
+        self.set_CFL(1.0)
+        self.set_use_kinematic_viscosity(False)
+
+        # Because gravity is treated within the flux function, 
+        # we remove it from the forcing terms.
+        self.forcing_terms.remove(gravity)
+        print 'Using shallow_water_balanced2 solver in /balanced_basic/'
+
+
+
+    #-----------------
+    # 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
+
+        # Shortcuts
+        Stage = domain.quantities['stage']
+        Xmom = domain.quantities['xmomentum']
+        Ymom = domain.quantities['ymomentum']
+        Bed = domain.quantities['elevation']
+
+        timestep = float(sys.maxint)
+
+        flux_timestep = compute_fluxes_ext(timestep,
+                                           domain.epsilon,
+                                           domain.H0,
+                                           domain.g,
+                                           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,
+                                           Stage.boundary_values,
+                                           Xmom.boundary_values,
+                                           Ymom.boundary_values,
+                                           Stage.explicit_update,
+                                           Xmom.explicit_update,
+                                           Ymom.explicit_update,
+                                           domain.already_computed_flux,
+                                           domain.max_speed,
+                                           int(domain.optimise_dry_cells),
+                                           Stage.centroid_values,
+                                           Bed.centroid_values, 
+                                           Bed.vertex_values)
+
+        #import pdb
+        #pdb.set_trace()
+
+        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
+        zc = domain.quantities['elevation'].centroid_values
+        zv = domain.quantities['elevation'].vertex_values
+        xmomc = domain.quantities['xmomentum'].centroid_values
+        ymomc = domain.quantities['ymomentum'].centroid_values
+
+        protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
+                domain.epsilon, wc, zc,zv, xmomc, ymomc)
+
+    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 """
+
+                #Shortcuts
+        #W  = self.quantities['stage']
+        #Z  = self.quantities['elevation']
+
+        #Arrays
+        #w_C   = W.centroid_values
+        #z_C   = Z.centroid_values
+
+        #num_min = num.min(w_C-z_C)
+        #if  num_min < 0.0:
+        #    print 'num.min(w_C-z_C)', num_min
+
+        #print 'using swb2 distribute_to_vertices_and_edges'
+        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_to_vertices_and_edges(self):
+    #    """Distribution from centroids to edges specific to the SWW eqn.
+
+    #    In addition, all conserved quantities get distributed as per either a
+    #    constant (order==1) or a piecewise linear function (order==2).
+    #    
+
+    #    Precondition:
+    #    All conserved quantities defined at centroids and bed elevation defined at
+    #    edges.
+    #    
+    #    Postcondition
+    #    Evolved quantities defined at vertices and edges
+    #    """
+    #    from swb2_domain import protect_against_infinitesimal_and_negative_heights
+
+    #    if self._order_ == 1:
+    #        for name in [ 'stage', 'xmomentum', 'ymomentum' ]:
+    #            Q = self.quantities[name]
+    #            Q.extrapolate_first_order()
+    #    elif self._order_ == 2:
+    #        for name in [ 'stage', 'xmomentum', 'ymomentum' ]:
+    #            Q = self.quantities[name]
+    #            #Q.extrapolate_second_order_and_limit_by_edge()
+    #            Q.extrapolate_second_order_and_limit_by_vertex()
+
+    #        #self.extrapolate_second_order_sw()
+    #    else:
+    #        raise Exception('Unknown order: %s' % str(self._order_))
+
+
+    #    # Update other quantities
+    #    protect_against_infinitesimal_and_negative_heights(self)
+
+    #    
+    #    # Compute edge values by interpolation
+    #    for name in ['stage', 'xmomentum', 'ymomentum']:
+    #        Q = self.quantities[name]
+    #        Q.interpolate_from_vertices_to_edges()
+
+
+    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
+        """
+        #from swb2_domain import protect_against_infinitesimal_and_negative_heights
+
+        #print 'using swb2 distribute_using_vertex_limiter'
+        # 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
+        domain.protect_against_infinitesimal_and_negative_heights()
+
+        #print 'using swb2 distribute_using_vertex_limiter'
+        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()
+            else:
+                raise Exception('Unknown order')
+
+        #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 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
+        """
+
+        # 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))
Index: /trunk/anuga_work/development/gareth/balanced_basic/swb2_domain_ext.c
===================================================================
--- /trunk/anuga_work/development/gareth/balanced_basic/swb2_domain_ext.c	(revision 8302)
+++ /trunk/anuga_work/development/gareth/balanced_basic/swb2_domain_ext.c	(revision 8302)
@@ -0,0 +1,1083 @@
+// 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 = 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;
+}
+
+// Innermost flux function (using stage w=z+h)
+int _flux_function_central(double *q_left, double *q_right,
+                           double z_left, double z_right,
+                           double n1, double n2,
+                           double epsilon, 
+                           double h0,
+                           double limiting_threshold, 
+                           double g,
+                           double *edgeflux, double *max_speed, 
+                           double *pressure_flux) 
+{
+
+  /*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, h_left, uh_left, vh_left, u_left;
+  double w_right, h_right, uh_right, vh_right, u_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double denom, inverse_denominator, z;
+  double uint, t1, t2, t3;
+  // 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 = 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 = w_left - z;
+  uh_left = q_left_rotated[1];
+  u_left = _compute_speed(&uh_left, &h_left, 
+			  epsilon, h0, limiting_threshold);
+
+  w_right = q_right_rotated[0];
+  h_right = w_right - z;
+  uh_right = q_right_rotated[1];
+  u_right = _compute_speed(&uh_right, &h_right, 
+			   epsilon, h0, limiting_threshold);
+
+  // 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;
+  }
+
+  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
+  if (s_min > 0.0)
+  {
+    s_min = 0.0;
+  }
+  
+  // Flux formulas
+  flux_left[0] = u_left*h_left;
+  flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
+  flux_left[2] = u_left*vh_left;
+
+  flux_right[0] = u_right*h_right;
+  flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
+  flux_right[2] = u_right*vh_right;
+
+  // Flux computation
+  denom = s_max - s_min;
+  if (denom < epsilon) 
+  { // FIXME (Ole): Try using h0 here
+    memset(edgeflux, 0, 3*sizeof(double));
+    //for(i=0;i<3;i++){
+    //    edgeflux[i] = _minmod(flux_left[i], flux_right[i]);
+    //}
+    *max_speed = 0.0;
+  } 
+  else 
+  {
+    inverse_denominator = 1.0/denom;
+    for (i = 0; i < 3; i++) 
+    {
+      // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
+      // Physics 2:141-163
+      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
+      //t1 = (q_right_rotated[i] - uint);
+      //t2 = (-q_left_rotated[i] + uint);
+      //t3 = _minmod(t1,t2);
+      t3 = 0.0;
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
+      //t1 = s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
+      //if(fabs(edgeflux[i])>fabs(t1)){ 
+      //      edgeflux[i]+=t1;
+      //}else{
+      //      edgeflux[i]+=fabs(edgeflux[i])*copysign(1.0,t1) ;
+      //}
+
+      edgeflux[i] *= inverse_denominator;
+    }
+
+    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
+    //if(edgeflux[0]<0.0){
+    //    //edgeflux[2] = edgeflux[0]*0.5*(vh_right/h_right + vh_left/h_left);
+    //    edgeflux[2] = edgeflux[0]*vh_right/h_right;
+    //}else{
+    //    edgeflux[2] = edgeflux[0]*vh_left/h_left;
+    //}
+
+    // Maximal wavespeed
+    *max_speed = max(fabs(s_max), fabs(s_min));
+
+    // 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,
+        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* stage_boundary_values,
+        double* xmom_boundary_values,
+        double* ymom_boundary_values,
+        double* stage_explicit_update,
+        double* xmom_explicit_update,
+        double* ymom_explicit_update,
+        long* already_computed_flux,
+        double* max_speed_array,
+        int optimise_dry_cells, 
+        double* stage_centroid_values,
+        double* bed_centroid_values,
+        double* bed_vertex_values) {
+    // Local variables
+    double max_speed, length, inv_area, zl, zr;
+    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
+
+    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
+    // threshold for performance reasons.
+    // See ANUGA manual under flux limiting
+    int k, i, m, n,j;
+    int ki, nm = 0, ki2; // 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;
+    int neighbours_wet[3];//Work array
+    int useint;
+    double usedouble, scale_factor_shallow, bedtop, bedbot, pressure_flux;
+    //double depthtemp, max_bed_edge_value; 
+    static long call = 1; // Static local variable flagging already computed flux
+
+    // Start computation
+    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));
+
+    // 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
+
+            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];
+
+            // Get right hand side values either from neighbouring triangle
+            // or from boundary array (Quantities at neighbour on nearest face).
+            n = neighbours[ki];
+            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
+            }
+            else {
+                // Neighbour is a real triangle
+                m = neighbour_edges[ki];
+                nm = n * 3 + m; // Linear index (triangle n, edge m)
+
+                qr[0] = stage_edge_values[nm];
+                qr[1] = xmom_edge_values[nm];
+                qr[2] = ymom_edge_values[nm];
+                zr = bed_edge_values[nm];
+            }
+
+            // Now we have values for this edge - both from left and right side.
+
+            if (optimise_dry_cells) {
+                // Check if flux calculation is necessary across this edge
+                // This check will exclude dry cells.
+                // This will also optimise cases where zl != zr as
+                // long as both are dry
+
+                if ((ql[0] - zl) < epsilon &&
+                        (qr[0] - zr) < epsilon) {
+                    // Cell boundary is dry
+
+                    already_computed_flux[ki] = call; // #k Done
+                    if (n >= 0) {
+                        already_computed_flux[nm] = call; // #n Done
+                    }
+
+                    max_speed = 0.0;
+                    continue;
+                }
+            }
+
+
+            if (fabs(zl-zr)>1.0e-10) {
+                report_python_error(AT,"Discontinuous Elevation");
+                return 0.0;
+            }
+            
+            // If both centroids are dry, then the cells should not exchange flux 
+            // NOTE: This also means no bedslope term -- which is arguably
+            // appropriate, since the depth is zero at the cell centroid
+            if( (stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
+                if(n>0){
+                    if(stage_centroid_values[n]-bed_centroid_values[n]<=0.0){
+                        continue;
+                    }
+                }else{
+                    // Treat boundaries -- if the boundary stage < zl, there is
+                    // no flux
+                    if(qr[0]< zl){
+                        continue;
+                    }    
+                }
+            }  
+           
+            //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
+            // unless the local centroid value is smaller 
+            if(n>=0){
+                if((stage_centroid_values[k]-bed_centroid_values[k]<= 0.0)){
+                          ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+                }
+                if(stage_centroid_values[n]-bed_centroid_values[n]<= 0.0){
+                        qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
+                }
+            }else{
+                // Treat the boundary case
+                if((stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
+                          ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+                }
+
+            }
+            
+            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
+            ki2 = 2 * ki; //k*6 + i*2
+
+            // Edge flux computation (triangle k, edge i)
+            _flux_function_central(ql, qr, zl, zr,
+                    normals[ki2], normals[ki2 + 1],
+                    epsilon, h0, limiting_threshold, g,
+                    edgeflux, &max_speed, &pressure_flux);
+            //_flux_function_fraccaro(ql, qr, zl, zr,
+            //        normals[ki2], normals[ki2 + 1],
+            //        epsilon, h0, limiting_threshold, g,
+            //        edgeflux, &max_speed);
+   
+            //Find the bed vertex values associated with the edge 
+            //if(i==0){
+            //    bedtop = max(bed_vertex_values[3*k+1], bed_vertex_values[3*k+2]);
+            //    bedbot = min(bed_vertex_values[3*k+1], bed_vertex_values[3*k+2]);
+            //}else{
+            //    if(i==1){
+            //    bedtop = max(bed_vertex_values[3*k+0], bed_vertex_values[3*k+2]);
+            //    bedbot = min(bed_vertex_values[3*k+0], bed_vertex_values[3*k+2]);
+            //    }else{
+            //    // i==2
+            //    bedtop = max(bed_vertex_values[3*k+1], bed_vertex_values[3*k+0]);
+            //    bedbot = min(bed_vertex_values[3*k+1], bed_vertex_values[3*k+0]);
+            //    }
+            //
+            //}       
+
+            //if(fabs(bedtop-bedbot)<1.0e-10){
+            //    scale_factor_shallow = 1.0;
+            //}else{
+            //    scale_factor_shallow = min( (bedtop - bedbot)/(min(stage_centroid_values[k], bedtop) - bedbot),10.0);
+            //}
+            
+            // Multiply edgeflux by edgelength
+            length = edgelengths[ki];
+            //length = edgelengths[ki]*scale_factor_shallow;
+            edgeflux[0] *= length;
+            edgeflux[1] *= length;
+            edgeflux[2] *= length;
+           
+            // To correctly implement reflective boundary conditions would be something like this? 
+            //if(n<0){
+            //   _rotate(edgeflux,normals[ki2],normals[ki2+1]);
+            //   edgeflux[1] = 0.0;
+            //   _rotate(edgeflux,-normals[ki2+1],normals[ki2]);
+            //}
+
+            // Update triangle k with flux from edge i
+            stage_explicit_update[k] -= edgeflux[0];
+            xmom_explicit_update[k] -= edgeflux[1];
+            ymom_explicit_update[k] -= edgeflux[2];
+
+            // Calculate the bed slope term, using the 'divergence form' from 
+            // Alessandro Valiani and Lorenzo Begnudelli, Divergence Form for Bed Slope Source Term in Shallow Water Equations
+            // Journal of Hydraulic Engineering, 2006, 132:652-665
+            if((stage_centroid_values[k]>bed_centroid_values[k]+0.0)){
+                bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
+                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
+                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
+            }else{
+                // Treat nearly dry cells 
+                //bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
+                bedslope_work = -pressure_flux*length;
+                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
+                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
+            }
+                // The triangle is (partially?) dry
+                // Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
+                // FIXME: Note that we keep recomputing this -- could improve efficiency 
+                //for(j=0;j<3;j++){ 
+                //    // For each edge, check if each neighbouring triangle is wet
+                //    useint = neighbours[3*k+j];
+                //    if(useint>0){
+                //    // neighbour is not a boundary condition
+                //        if(stage_centroid_values[useint] > bed_centroid_values[useint]){
+                //            //wet
+                //            neighbours_wet[j]=1;
+                //            useint = useint*3 + neighbour_edges[3*k+j]; // Index of neighbouring edge
+                //            stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
+                //        }else{
+                //            //dry
+                //            neighbours_wet[j]=0;
+                //            stage_edges[j]=0.0; // Arbitrary
+                //        }
+                //    }else{ 
+                //        // neighbour is a boundary condition, ASSUME IT IS WET
+                //        useint = -useint -1;
+                //        neighbours_wet[j]=1;
+                //        stage_edges[j] = stage_boundary_values[useint];
+                //        }         
+                //}
+                ////calculate average neighbouring edge stage, only considering wet triangles
+                //if(max(neighbours_wet[0], max(neighbours_wet[1],neighbours_wet[2]))>0){
+                //    //usedouble == average edge stage
+                //    usedouble = stage_edges[0]*neighbours_wet[0] + stage_edges[1]*neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
+                //    usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
+                //}else{
+                //    usedouble = zl;
+                //}
+                //usedouble=min(usedouble,stage_centroid_values[k]);
+                //xmom_explicit_update[k] -= -0.5*g*max(usedouble-zl,0.0)*(usedouble-zl)*normals[ki2]*length;
+                //ymom_explicit_update[k] -= -0.5*g*max(usedouble-zl,0.0)*(usedouble-zl)*normals[ki2+1]*length;
+                ////xmom_explicit_update[k] -= -0.5*g*abs(usedouble-zl)*(usedouble-zl)*normals[ki2]*length;
+                ////ymom_explicit_update[k] -= -0.5*g*abs(usedouble-zl)*(usedouble-zl)*normals[ki2+1]*length;
+
+            //}
+
+            already_computed_flux[ki] = call; // #k Done
+
+
+            // Update neighbour n with same flux but reversed sign
+            if (n >= 0) {
+                stage_explicit_update[n] += edgeflux[0];
+                xmom_explicit_update[n] += edgeflux[1];
+                ymom_explicit_update[n] += edgeflux[2];
+                //Add bed slope term here
+                if((stage_centroid_values[n]>bed_centroid_values[n]+0.0)){
+                    bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
+                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
+                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
+                }else{
+                    // Treat nearly dry cells
+                    //bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
+                    bedslope_work = -pressure_flux*length;
+                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
+                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
+                }
+                    // The triangle is (partially?) dry
+                    // Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
+                    // FIXME: Note that we keep recomputing this -- could improve efficiency
+                    //for(j=0;j<3;j++){
+                    //    // For each edge, check if each neighbouring triangle is wet
+                    //    useint = neighbours[3*n+j];
+                    //    if(useint>0){
+                    //    // neighbour is not a boundary condition
+                    //        if(stage_centroid_values[useint] > bed_centroid_values[useint]){
+                    //            //wet
+                    //            neighbours_wet[j]=1;
+                    //            useint = useint*3 + neighbour_edges[3*n+j]; // Index of neighbouring edge
+                    //            stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
+                    //        }else{
+                    //            //dry
+                    //            neighbours_wet[j]=0;
+                    //            stage_edges[j]=0.0; //Arbitrary
+                    //        }
+                    //    }else{ 
+                    //        // neighbour is a boundary condition, ASSUME IT IS WET
+                    //        useint = -useint -1;
+                    //        neighbours_wet[j]=1;
+                    //        stage_edges[j] = stage_boundary_values[useint];
+                    //        }         
+                    //}
+                    ////calculate average neighbouring edge stage, only considering wet triangles
+                    //if(max(neighbours_wet[0],max(neighbours_wet[1],neighbours_wet[2]))>0){
+                    //    //average edge stage
+                    //    usedouble = stage_edges[0]*neighbours_wet[0] + stage_edges[1]*neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
+                    //    usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
+                    //}else{
+                    //    usedouble = zr;
+                    //}
+                    //usedouble=min(usedouble,stage_centroid_values[n]);
+                    //xmom_explicit_update[n] += -0.5*g*max(usedouble-zr,0.0)*(usedouble-zr)*normals[ki2]*length;
+                    //ymom_explicit_update[n] += -0.5*g*max(usedouble-zr,0.0)*(usedouble-zr)*normals[ki2+1]*length;
+                    ////xmom_explicit_update[n] += -0.5*g*abs(usedouble-zr)*(usedouble-zr)*normals[ki2]*length;
+                    ////ymom_explicit_update[n] += -0.5*g*abs(usedouble-zr)*(usedouble-zr)*normals[ki2+1]*length;
+
+               //}
+
+                already_computed_flux[nm] = call; // #n Done
+            }
+
+            // Update timestep based on edge i and possibly neighbour n
+            if (tri_full_flag[k] == 1) {
+                if (max_speed > epsilon) {
+                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
+
+                    // CFL for triangle k
+                    timestep = min(timestep, radii[k] / max_speed);
+
+                    if (n >= 0) {
+                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
+                        timestep = min(timestep, radii[n] / max_speed);
+                    }
+
+                    // Ted Rigby's suggested less conservative version
+                    //if(n>=0){
+                    //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
+                    //}else{
+                    //  timestep = min(timestep, radii[k]/max_speed);
+                    //}
+                }
+            }
+
+        } // End edge i (and neighbour n)
+
+
+        // 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;
+
+
+        // Keep track of maximal speeds
+        max_speed_array[k] = max_speed;
+
+    } // End triangle k
+
+
+    // Try to adjust stage_explicit_update to account for the negative stages,
+    // using an approach assumes a very small discharge is enough to allow 
+    // the depth to increase from negative to zero
+    //for(k=0;k<number_of_elements;k++){
+        //depthtemp = stage_centroid_values[k]-bed_centroid_values[k];
+    //    depthtemp = stage_centroid_values[k]-max(stage_edge_values[3*k+1], max(stage_edge_values[3*k+2], stage_edge_values[3*k]));
+    //   if(depthtemp< 0.0){
+    //       depthtemp = -depthtemp; // stage defecit
+    //       usedouble = 2.0; //1.0e+2; //Large constant
+    //       if(stage_explicit_update[k]*timestep*usedouble > depthtemp){
+    //            stage_explicit_update[k] += (depthtemp - depthtemp/usedouble)/timestep;
+    //        }else{
+    //            stage_explicit_update[k] *=usedouble ;
+    //        } 
+    //    }
+    //}
+
+    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* zc,
+         double* zv,
+         double* xmomc,
+         double* ymomc) {
+
+  int k;
+  double hc, wmin;
+  double u, v, reduced_speed;
+
+
+  // Protect against initesimal and negative heights  
+  if (maximum_allowed_speed < epsilon) {
+    for (k=0; k<N; k++) {
+      hc = wc[k] - zc[k];
+
+      if (hc < minimum_allowed_height) {
+        
+        // Set momentum to zero and ensure h is non negative
+        xmomc[k] = 0.0;
+        ymomc[k] = 0.0;
+        if (hc <= 0.0){
+             // Definine the minimum possible stage value as the minimum bed edge value on triangle k.
+             wmin = 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])));
+             wc[k] = max(wc[k], wmin);
+        }
+      }
+    }
+  } else {
+   
+    // Protect against initesimal and negative heights
+    for (k=0; k<N; k++) {
+      hc = wc[k] - zc[k];
+      
+      if (hc < minimum_allowed_height) {
+
+        //New code: Adjust momentum to guarantee speeds are physical
+        //          ensure h is non negative
+              
+        if (hc <= 0.0) {
+            wc[k] = zc[k];
+        xmomc[k] = 0.0;
+        ymomc[k] = 0.0;
+        } else {
+          //Reduce excessive speeds derived from division by small hc
+	  //FIXME (Ole): This may be unnecessary with new slope limiters 
+	  //in effect.
+          
+          u = xmomc[k]/hc;
+	  if (fabs(u) > maximum_allowed_speed) {
+	    reduced_speed = maximum_allowed_speed * u/fabs(u);
+	    //printf("Speed (u) has been reduced from %.3f to %.3f\n",
+	    //   u, reduced_speed);
+	    xmomc[k] = reduced_speed * hc;
+	  }
+	  
+          v = ymomc[k]/hc;
+	  if (fabs(v) > maximum_allowed_speed) {
+	    reduced_speed = maximum_allowed_speed * v/fabs(v);
+	    //printf("Speed (v) has been reduced from %.3f to %.3f\n",
+	    //   v, reduced_speed);
+	    ymomc[k] = reduced_speed * hc;
+	  }
+        }
+      }
+    }
+  }
+  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 *neighbours, *neighbour_edges,
+    *normals, *edgelengths, *radii, *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, //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,
+    *bed_centroid_values,
+    *bed_vertex_values;
+    
+  double timestep, epsilon, H0, g;
+  int optimise_dry_cells;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOiOOO",
+            &timestep,
+            &epsilon,
+            &H0,
+            &g,
+            &neighbours,
+            &neighbour_edges,
+            &normals,
+            &edgelengths, &radii, &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,
+            &max_speed_array,
+            &optimise_dry_cells,
+            &stage_centroid_values,
+            &bed_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(stage_boundary_values);
+  CHECK_C_CONTIG(xmom_boundary_values);
+  CHECK_C_CONTIG(ymom_boundary_values);
+  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(bed_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,
+                     (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*) stage_boundary_values -> data,
+                     (double*) xmom_boundary_values -> data,
+                     (double*) ymom_boundary_values -> 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, 
+                     (double*) stage_centroid_values -> data, 
+                     (double*) bed_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;
+
+  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;
+  }
+
+  
+  h0 = H0*H0; // This ensures a good balance when h approaches H0.
+              // But evidence suggests that h0 can be as little as
+              // epsilon!
+	      
+  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 
+  _flux_function_central((double*) ql -> data, 
+			 (double*) qr -> data, 
+			 zl, 
+			 zr,                         
+			 ((double*) normal -> data)[0],
+			 ((double*) normal -> data)[1],          
+			 epsilon, h0, limiting_threshold,
+			 g,
+			 (double*) edgeflux -> data, 
+			 &max_speed,
+             &pressure_flux);
+  
+  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("");
+}
+
+
+//========================================================================
+// 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
+  *zc,            //Elevation at centroids
+  *zv,            //Elevation at vertices
+  *xmomc,         //Momentums at centroids
+  *ymomc;
+
+
+  int N;
+  double minimum_allowed_height, maximum_allowed_speed, epsilon;
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dddOOOOO",
+            &minimum_allowed_height,
+            &maximum_allowed_speed,
+            &epsilon,
+            &wc, &zc, &zv, &xmomc, &ymomc)) {
+    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*) zc -> data,
+       (double*) zv -> data,
+       (double*) xmomc -> data,
+       (double*) ymomc -> 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.
+   */
+  {"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"},  
+  {"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
+}
