source: trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

// 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 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, double hc, 
                           double hc_n, double speed_max_last, 
                           double smooth) 
{

  /*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, 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 = 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.);
  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.);
  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];

      // Standard smoothing term
      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);

      // Standard smoothing term, scaled
      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*smooth;

      // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
      // Physics 2:141-163
      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
      //t1 = (q_right_rotated[i] - uint);
      //t2 = (-q_left_rotated[i] + uint);
      //t3 = _minmod(t1,t2);
      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3);

      // test this
      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*pow(min(hc/(hc_n+epsilon), hc_n/(hc+epsilon)),1.0);
    
      // test this 
      //tmp=min(fabs(q_right_rotated[i])/(fabs(q_left_rotated[i])+epsilon) , fabs(q_left_rotated[i])/(fabs(q_right_rotated[i])+epsilon));
      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3)*pow(tmp,1.0);

      //GD HACK: Here, we supress the 'smoothing' part of the flux -- like scaling diffusion by local wave speed
      // This is a bad idea! E.g. oscillations in landslide wave runup case
      // However, it does supress the growth of some wet-dry instabilities (e.g. PNG).
      //if(i!=2){
      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*
      //               (min(*max_speed/(max(speed_max_last,1.0e-100)),1.0));//(min(hc/h_left,min(hc_n/h_right,1.0)));
      //}

      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* 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* 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 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 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, *count_wet_neighbours, *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 minimum bed edge value on each triangle
    speed_max_last=0.0;
    for (k = 0; k < number_of_elements; k++){
        
        count_wet_edges[k] = 0.;
        count_wet_neighbours[k] = 0.;
        for (i =0; i<3; i++){
            ki=3*k+i;
            if(stage_edge_values[ki] > bed_edge_values[ki]+epsilon) count_wet_edges[k]+=1.0;
            
            n=neighbours[ki];
            if(n<0 || stage_centroid_values[n] > bed_centroid_values[n]+h0+epsilon) count_wet_neighbours[k]+=1.0;
        }

        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 = max(stage_centroid_values[k] - bed_centroid_values[k],0.0);

            // 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;
            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
                hc_n = max(stage_centroid_values[n] - bed_centroid_values[n],0.0);
                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];
            }
            
            if (fabs(zl-zr)>1.0e-10) {
                report_python_error(AT,"Discontinuous Elevation");
                return 0.0;
            }
           

            smooth=1.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, zl, zr,
                    normals[ki2], normals[ki2 + 1],
                    epsilon, h0, limiting_threshold, g,
                    edgeflux, &max_speed, &pressure_flux, hc, hc_n, 
                    speed_max_last, smooth);

           
            // 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.

            if(ql[0]>zl){
               tmp=1.0;
            }else{
               tmp=1.;
            }
            if(hc> h0){
            //if(stage_centroid_values[k] > count_wet_edges[k]){
              bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl))*tmp + pressure_flux);
              //bedslope_work = length*(g*(hc*stage_centroid_values[k]*0.-
              //                           0.5*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)) 
              //                        + 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]*tmp-0.0*max(ql[0]-zl,0.)*(ql[0]-zl))*tmp + 0.0*pressure_flux);
              //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
              //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
              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];

                if(qr[0]>zr){
                   tmp=1.0;
                }else{
                   tmp=1.;
                }
                if(hc_n> h0){
                //if(stage_centroid_values[n] > count_wet_edges[n]){
                  bedslope_work = length*(g*(hc_n*qr[0]*tmp-0.5*max(qr[0]-zr,0.)*(qr[0]-zr))*tmp + pressure_flux);
                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]*0.-0.5*max(stage_centroid_values[n]-zr,0.)*
                  //                                (stage_centroid_values[n]-zr)) + pressure_flux);
                  //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]*tmp-0.0*max(qr[0]-zr,0.)*(qr[0]-zr))*tmp + 0.0*pressure_flux);
                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
                  //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
                  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* edge_coordinates,
                                 double* stage_edge_values,
                                 double* xmom_edge_values,
                                 double* ymom_edge_values,
                                 double* elevation_edge_values,
                                 double* stage_vertex_values,
                                 double* xmom_vertex_values,
                                 double* ymom_vertex_values,
                                 double* elevation_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, *stage_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));
  stage_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]));
            
          }

      }

  //
  // 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(count_wet_neighbours[k]<=3){
    if(1==1){
        bedmax = max(elevation_vertex_values[3*k],
                         max(elevation_vertex_values[3*k+1], 
                             elevation_vertex_values[3*k+2]));
        if((cell_wetness_scale[k]==0. ) ||
            ((stage_centroid_values[k] - elevation_centroid_values[k] < 0.05*(bedmax-elevation_centroid_values[k])) &
             count_wet_neighbours[k]<=3) ){
            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];
          
          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))
          {

              if(0==1){
                //// Friction slope type extrapolation -- emulating steady flow
                //// d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth**(4/3)
                //hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
                //if(hc>1.0e-06){
                //  tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
                //             ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
                //  tmp /=pow(hc,10.0/3.0);
                //}else{
                //  tmp=0.0;
                //}
                //a = -xmom_centroid_values[k]*tmp;
                //b = -ymom_centroid_values[k]*tmp;
            
                //// Limit gradient to be between the bed slope, and zero
                //if(a*a_bs<0.) a=0.;
                //if(fabs(a)>fabs(a_bs)) a=a_bs;
                //if(b*b_bs<0.) b=0.;
                //if(fabs(b)>fabs(b_bs)) b=b_bs;


                //dqv[0] = a*dxv0 + b*dyv0;
                //dqv[1] = a*dxv1 + b*dyv1;
                //dqv[2] = a*dxv2 + b*dyv2;

                //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];
              }else{
                // Constant stage extrapolation
                //stage_edge_values[k3]   = stage_centroid_values[k];
                //stage_edge_values[k3+1] = stage_centroid_values[k];
                //stage_edge_values[k3+2] = stage_centroid_values[k];
              }
            //}else{
            //}else{
              // Dry cell with a single 'wet' neighbour

              // Compute indices of neighbouring wet cell / edge
              //ktmp = k0*(1-(k0==k)) + k1*(1-(k1==k)) + k2*(1-(k2==k));
              //ii = 0*(1-(k0==k)) + 1*(1-(k1==k)) + 2*(1-(k2==k));
              //if((ii<0)|(ii>2)){
              //  printf("Invalid ii value \n");
              //}
            
              //ktmp = 3*ktmp+neighbour_edges[3*k+ii];
              //stage_edge_values[k3]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
              //stage_edge_values[k3+1]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
              //stage_edge_values[k3+2]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);

            //}   
              if( (k==k0) & (k==k1) & (k==k2)){
              //    // Isolated wet cell  -- use constant depth extrapolation for stage
              //    stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
              //    stage_edge_values[k3+1] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1];
              //    stage_edge_values[k3+2] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2];
              //}
              //    // Isolated wet cell -- constant depth extrapolation, but don't let water flow uphill
              //    //stage_edge_values[k3]   = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3]);
              //    //stage_edge_values[k3+1] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1]);
              //    //stage_edge_values[k3+2] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2]);
              //    // Isolated wet cell -- constant stage extrapolation 
                  stage_edge_values[k3]   = stage_centroid_values[k];
                  stage_edge_values[k3+1] = stage_centroid_values[k];
                  stage_edge_values[k3+2] = stage_centroid_values[k];
              }else{
              //    // Cell with one wet neighbour. 
              //    // Find which neighbour is wet [if k!=k0, then k0 is wet]
                  if(k!=k0){
                    ktmp=k0;
                    ii=0;
                    xtmp = x0;
                    ytmp = y0;
                  }else if(k!=k1){
                    ktmp=k1;
                    ii=1;
                    xtmp = x1;
                    ytmp = y1;
                  }else if(k!=k2){
                    ktmp=k2;
                    ii=2;
                    xtmp = x2;
                    ytmp = y2;
                  }else{
                      printf("ERROR in extrapolation routine: Should have had one ki == k \n");
                      report_python_error(AT, "Negative triangle area");
                      return -1;
                  }

              //    //if(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
              //    if(k_wetdry==0){
              //    //    // Cell is wet
              //    //    //if(count_wet_neighbours[ktmp]>0){
              //    //    //    stage_edge_values[k3]= stage_centroid_values[k];
              //    //    //    stage_edge_values[k3+1]= stage_centroid_values[k];
              //    //    //    stage_edge_values[k3+2]= stage_centroid_values[k];
              //    //    //}else{
              //    //        // Constant depth extrapolation
                    if(stage_centroid_values[k]<stage_centroid_values[ktmp]){
                      // Constant stage extrapolation reduces the issue of having the bed-slope term playing up on nearly dry cells
                      // stage_edge_values[k3]= stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3];
                      // stage_edge_values[k3+1]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1];
                      // stage_edge_values[k3+2]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2];
                      stage_edge_values[k3] = stage_centroid_values[k];
                      stage_edge_values[k3+1] = stage_centroid_values[k];
                      stage_edge_values[k3+2] = stage_centroid_values[k];
                    }else{
              //           // stage_edge_values[k3] = stage_centroid_values[ktmp];
              //           // stage_edge_values[k3+1] = stage_centroid_values[ktmp];
              //           // stage_edge_values[k3+2] = stage_centroid_values[ktmp];
                     //stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
                     //stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
                     //stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
                     ktmp = 3*ktmp+neighbour_edges[k3+ii];
                     stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_edge_values[ktmp]);
                     stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_edge_values[ktmp]);
                     stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_edge_values[ktmp]);
                    }

              //    }else{
              //        
              //    //    // This cell is 'dry'.
              //    //    // Extrapolate using the neighbouring wet cell if appropriate
              //    //    // This needs to preserve a wet-dry interface
              //        //ktmp = 3*ktmp+neighbour_edges[3*k+ii];
              //        //stage_edge_values[k3]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
              //        //stage_edge_values[k3+1]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
              //        //stage_edge_values[k3+2]= stage_edge_values[ktmp];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
              //    //    stage_edge_values[k3]= stage_centroid_values[ktmp] ;
              //    //    stage_edge_values[k3+1]= stage_centroid_values[ktmp];
              //    //    stage_edge_values[k3+2]= stage_centroid_values[ktmp];
              //        stage_edge_values[k3]= stage_centroid_values[k] ;
              //        stage_edge_values[k3+1]= stage_centroid_values[k];
              //        stage_edge_values[k3+2]= stage_centroid_values[k];
              //    }
              //
              }
              // First order momentum / velocity extrapolation
             //if(stagemin>=bedmax){
             //   stage_edge_values[k3]    = stage_centroid_values[k];
             //   stage_edge_values[k3+1]  = stage_centroid_values[k];
             //   stage_edge_values[k3+2]  = stage_centroid_values[k];
             //}else{
             //   stage_edge_values[k3]= stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3];
             //   stage_edge_values[k3+1]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1];
             //   stage_edge_values[k3+2]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2];
             //}
              xmom_edge_values[k3]    = xmom_centroid_values[k];
              xmom_edge_values[k3+1]  = xmom_centroid_values[k];
              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(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height*0.1), 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];
         
          //if( (area2>0.0) ){ 

          inv_area2 = 1.0/area2;
          //if(count_wet_neighbours[k] > 1){
          //if(1==1){
          //    // 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;
          //}else{
          ////    //// Friction slope type extrapolation -- emulating steady flow
          ////    //// d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth**(4/3)
          //    hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
          //    if(hc>1.0e-06){
          //      tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
          //               ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
          //      tmp /=max(pow(hc,10.0/3.0), 1.0e-30);
          //    }else{
          //      tmp=0.0;
          //    }
          //    a = -xmom_centroid_values[k]*tmp;
          //    b = -ymom_centroid_values[k]*tmp;
          ////    // Limit gradient to be between the bed slope, and zero
          ////    if(a*a_bs<0.) a=0.;
          ////    if(fabs(a)>fabs(a_bs)) a=a_bs;
          ////    if(b*b_bs<0.) b=0.;
          ////    if(fabs(b)>fabs(b_bs)) b=b_bs;
          //}
          // 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
          // This is more complex for stage than other quantities
          //if((count_wet_neighbours[k]>0)&&(cell_wetness_scale[k]=1.0)){//&&(k_wetdry==0)){
          //if(stage_centroid_values[k] > max_elevation_edgevalue[k]){
          //if( (area2>0) ){ 
          //if(count_wet_neighbours[k]>0){
          //if(min_elevation_edgevalue[k]<=stagemin ){
          //if(bedmax<=stagemin ){
             limit_gradient(dqv, qmin, qmax, beta_tmp);
          //}else{
          //  dqv[0]=-elevation_centroid_values[k] + elevation_edge_values[k3+0];
          //  dqv[1]=-elevation_centroid_values[k] + elevation_edge_values[k3+1];
          //  dqv[2]=-elevation_centroid_values[k] + elevation_edge_values[k3+2];
          //}
          //printf("%lf \n ", beta_tmp);
          //}
          //}
          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];
          // IDEA: extrapolate assuming slope is that of steady flow??
          // GD HACK - FIXME
          // Note that 'near-dry' cells which don't exchange mass with
          // neighbouring cells have to revert to having a constant stage
          // extrapolation, or else they will have a residual bed slope term.
          // This helps! -- But probably makes rainfall-runoff problems worse.
          //ii = (stage_edge_values[k3+0] > elevation_edge_values[k3+0])+ 
          //     (stage_edge_values[k3+1] > elevation_edge_values[k3+1])+ 
          //     (stage_edge_values[k3+2] > elevation_edge_values[k3+2]);
          //if(ii<2){
          //    stage_edge_values[k3+0] = stage_centroid_values[k] ;
          //    stage_edge_values[k3+1] = stage_centroid_values[k] ;
          //    stage_edge_values[k3+2] = stage_centroid_values[k] ;
          //}
          //}else{
          // d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth**(4/3)
          //    hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
          //    tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
          //               ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
          //    tmp /=max(hc*hc*hc, 1.0e-09);
          //    a = -xmom_centroid_values[k]*tmp;
          //    b = -ymom_centroid_values[k]*tmp;
          //    dqv[0] = a*dxv0 + b*dyv0;
          //    dqv[1] = a*dxv1 + b*dyv1;
          //    dqv[2] = a*dxv2 + b*dyv2;
          //////if(count_wet_neighbours[k]==0){
          ////    //weight=max(cell_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
          ////    //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
          ////    //weight=min(weight, 1.0);
          ////    //weight==0;
          //      weight=0.;
          ////    // 'Shallow flow on steep slope'
              //limit_gradient(dqv, qmin, qmax, beta_tmp);
              //stage_edge_values[k3+0] =stage_centroid_values[k] + dqv[0];
              //stage_edge_values[k3+1] =stage_centroid_values[k] + dqv[1];
              //stage_edge_values[k3+2] =stage_centroid_values[k] + dqv[2];
          //
          ////    // There exists a neighbouring bed cell with elevation > minimum
          ////    // neighbouring stage value

          ////    // We have to be careful in this situation. Limiting the stage gradient can
          ////    // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
          ////    // Previously 'balance_deep_and_shallow' was used to deal with this issue
          //    for(i=0; i<3; i++){
          //        stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] + 
          //                                 (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
          //                                               +elevation_edge_values[k3+i]);
          //    }
          //}
           
          //-----------------------------------
          // 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);
          
          //for (i=0; i < 3; i++)
          //{
          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];
          //}
          
          //-----------------------------------
          // 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 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] = 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] = 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];
          }

      }

   
  } 

  free(xmom_centroid_store);
  free(ymom_centroid_store);
  free(stage_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,
    *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,
    *bed_vertex_values;
    
  double timestep, epsilon, H0, g;
  int optimise_dry_cells, timestep_order;
    
  // Convert Python arguments to C
  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOiiOOOOO",
            &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,
            &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,
            &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(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(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*) 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*) 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;

  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;            
  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,                         
                         ((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],
             smooth 
             );
  
  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,
    *edge_coordinates,
    *stage_edge_values,
    *xmom_edge_values,
    *ymom_edge_values,
    *elevation_edge_values,
    *stage_vertex_values,
    *xmom_vertex_values,
    *ymom_vertex_values,
    *elevation_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, "OOOOOOOOOOOOOOOOOOii",
            &domain,
            &surrogate_neighbours,
            &neighbour_edges,
            &number_of_boundaries,
            &centroid_coordinates,
            &stage_centroid_values,
            &xmom_centroid_values,
            &ymom_centroid_values,
            &elevation_centroid_values,
            &edge_coordinates,
            &stage_edge_values,
            &xmom_edge_values,
            &ymom_edge_values,
            &elevation_edge_values,
            &stage_vertex_values,
            &xmom_vertex_values,
            &ymom_vertex_values,
            &elevation_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(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(stage_vertex_values);
  CHECK_C_CONTIG(xmom_vertex_values);
  CHECK_C_CONTIG(ymom_vertex_values);
  CHECK_C_CONTIG(elevation_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*) edge_coordinates -> data,
                   (double*) stage_edge_values -> data,
                   (double*) xmom_edge_values -> data,
                   (double*) ymom_edge_values -> data,
                   (double*) elevation_edge_values -> data,
                   (double*) stage_vertex_values -> data,
                   (double*) xmom_vertex_values -> data,
                   (double*) ymom_vertex_values -> data,
                   (double*) elevation_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
}