1 0 GfsWave GfsBox GfsGEdge {} {
    Refine 4

    # Run for 48 hours
    Time { end = 48 }

    # Domain size is 3328 km
    PhysicalParams { L = 3328 }

    # Define some useful functions
    Global {
        /* gaussian distribution */
        static double gaussian (double f, double fmean, double fsigma) {
            return exp (-((f - fmean)*(f - fmean))/(fsigma*fsigma));
        }
        /* cos(theta)^n distribution */
        static double costheta (double theta, double thetam, double thetapower) {
            double a = cos (theta - thetam);
            return a > 0. ? pow (a, thetapower) : 0.;
        }
        /* Holland cyclone model */
        static double holland(double r, double Rmax, double Vmax) {
            if (r < Rmax/1e3) return 0.;
            return Vmax*pow(Rmax/r, 2)*exp(2.*(1. - Rmax/r));
        }
        /* Position of the center of the cyclone */
        double ut = 555./24.; /* km/h */
        static double xc (double t) {
            return 0.;
        }
        static double yc (double t) {
            return 1110. - ut*t;
        }
        /* Intensity of the cyclone as a function of time */
        static double vmax (double t) {
            return 50.*(t < 25. ? t/25. : 1.);
        }
        /* velocity components */
        static double ur (double x, double y, double t) {
            x -= xc (t);
            y -= yc (t);
            double r = sqrt (x*x + y*y);
            return holland (r, 100., vmax (t))*y/r;
        }
        static double vr (double x, double y, double t) {
            x -= xc (t);
            y -= yc (t);
            double r = sqrt (x*x + y*y);
            return - holland (r, 100., vmax (t))*x/r;
        }
    }

    # Use source terms from WaveWatch III
    GModule wavewatch

    Init { istep = 1 } {
        # Wind at 10 metres
        U10 = ur(x, y, t)
        V10 = vr(x, y, t)
    }

    # Adapt the mesh according to the error in significant wave height
    AdaptError { istep = 1 } { cmax = 0.1 minlevel = 4 maxlevel = LEVEL c = Hse } Hs
    # Adapt the mesh according to the error in the norm of the forcing wind field
    AdaptError { istep = 1 } { 
        cmax = 0.2 minlevel = 4 maxlevel = LEVEL c = Ve 
    } sqrt(U10*U10 + V10*V10)

    # Output time at every timestep
    OutputTime { istep = 1 } stderr
    # Output simulation size
    OutputBalance { istep = 1 } stderr
    # Output timing statistics every 100 timesteps
    OutputTiming { istep = 100 } stderr

    # Output significant wave height to file hs, every quarter of an hour
    OutputScalarStats { step = 0.25 } hs { v = Hs }
    # Output norm of wind velocity to file vr, every quarter of an hour
    OutputScalarStats { step = 0.25 } vr { v = sqrt(U10*U10 + V10*V10) }

    # Output simulation to standard output for visualisation with GfsView
    OutputSimulation { istep = 10 } stdout
    # Output simulation results every 4 hours
    OutputSimulation { step = 4 } sim-%g.gfs
    # Compress the files to save disk space
    EventScript { step = 4 } { gzip -f sim-*.gfs }
    # Create movies of significant wave height and level of refinement
    OutputPPM { istep = 1 } { ppm2mpeg > hs.mpg } { v = Hs maxlevel = 9 }
    OutputPPM { istep = 1 } { ppm2mpeg > level.mpg } { v = Level min = 4 max = LEVEL maxlevel = 9 }

    # Create figures at the end of the simulation
    EventScript { start = end } {
	for i in 12 24 36 48; do
	    echo "Save hs-$i.eps { format = EPS }" | gfsview-batch2D sim-$i.gfs.gz hs.gfv
	    echo "Save mesh-$i.eps { format = EPS }" | gfsview-batch2D sim-$i.gfs.gz mesh.gfv
	done
	echo "Save shifted.eps { format = EPS }" | gfsview-batch2D sim-48.gfs.gz shifted.gfv
	cat <<EOF | gnuplot
        set term postscript eps color lw 2 18
        set output 'hsmax.eps'
        set xlabel 'Time (hours)'
        set ylabel 'Amplitude (m or m/s)'
        set xtics 0,12,48
        set grid
        plot 'hs' u 3:11 w l t 'max(Hs)', 'vr' u 3:11 w l t 'max(|U10|)'
EOF
    }
} {
    # Garden Sprinkler Effect alleviation parameter
    alpha_s = ALPHA
}
GfsBox {}