Define h0 10.
Define a 3000.
Define tau 1e-3
Define B 5
Define G 9.81

1 0 GfsRiver GfsBox GfsGEdge {} {

    # Analytical solution, see Sampson, Easton, Singh, 2006
    Global {
	static gdouble Psi (double x, double t) {
	    double p = sqrt (8.*G*h0)/a;
	    double s = sqrt (p*p - tau*tau)/2.;
	    return a*a*B*B*exp (-tau*t)/(8.*G*G*h0)*(- s*tau*sin (2.*s*t) + 
		(tau*tau/4. - s*s)*cos (2.*s*t)) - B*B*exp(-tau*t)/(4.*G) -
	    exp (-tau*t/2.)/G*(B*s*cos (s*t) + tau*B/2.*sin (s*t))*x;
	}
    }

    PhysicalParams { L = 10000 }
    RefineSurface LEVEL (y - 10000.*(0.5 - 1./pow(2,LEVEL)))
    InitMask {} (y < 10000.*(0.5 - 1./pow(2,LEVEL)))
    Init {} {
	Zb = h0*(x/a)*(x/a)
	P = MAX (0., h0 + Psi (x, 0.) - Zb)
    }
    Init { istep = 1 } {
	Pt = MAX (0., h0 + Psi (x, t) - Zb)
    }
    # Better convergence rates are obtained at lower CFL
    # AdvectionParams { cfl = 0.1 }
    PhysicalParams { g = G }
    SourceCoriolis 0 tau
    Time { end = 6000 }
    OutputSimulation { start = 1500 } sim-LEVEL-%g.txt { format = text }
    OutputSimulation { start = end } end-LEVEL.txt { format = text }
    OutputScalarSum { istep = 10 } ke-LEVEL { v = (P > 0. ? U*U/P : 0.) }
    OutputScalarSum { step = 50 } vol-LEVEL { v = P }
    OutputScalarSum { step = 50 } U-LEVEL { v = U }
    OutputErrorNorm { istep = 1 } error-LEVEL { v = P } {
	s = Pt
	v = DP
    }
} {
    # this is necessary to obtain good convergence rates at high
    # resolutions
    dry = 1e-4
}
GfsBox {
    left = Boundary
    right = Boundary
}