     double uhat = (ur*hrs)/(hrs+hls) + (ul*hls)/(hrs + hls);
     double hhat = (hr + hl)/2;
 
-    double roe1 = vhat - sqrt(grav*hhat);
-    double roe3 = vhat + sqrt(grav*hhat);
-
-    double s1l = vl - sqrt(grav*hl);
-    double s3r = vr + sqrt(grav*hr);
+    double roe1 = vhat - sqrt(s_grav*hhat);
+    double roe3 = vhat + sqrt(s_grav*hhat);
+
+    double s1l = vl - sqrt(s_grav*hl);
+    double s3r = vr + sqrt(s_grav*hr);
 
     double s1 = min(roe1,s1l);
     double s3 = max(roe3,s3r);



        /* 1 wave speed of left state */
        double sl = ul - sqrt(grav*hl);    

        /* 2 wave speed of right state */
        double sr = ur + sqrt(grav*hr);    

        /* Roe speeds */
        double uhat = (sqrt(grav*hl)*ul + sqrt(grav*hr)*ur)/(sqrt(grav*hr) + sqrt(grav*hl));
        double chat = sqrt(grav*(hr + hl))/2;

        /* Compute wave speeds by comparing speeds above */
        s[0] = min(sl, uhat - chat);

        
        /* Flux differences */
        double fluxdiff[3];
        double phir = s_grav*hr*hr/2 + hur*hur/hr;
        double phil = s_grav*hl*hl/2 + hul*hul/hl;

        fluxdiff[0] = (hr*ur) - (hl*ul);
        fluxdiff[1] = phir - phil + s_grav*hbar*(br - bl);
        fluxdiff[2] = hr*ur*vr - hl*ul*vl;

        /* Wave strengths */
