# M2 tidal frequency. The period is 12h25 (44700 seconds).
Define M2F (2.*M_PI/44700.)

# M2 tidal elevation
Define M2(t) (A_amp*cos (M2F*t)+B_amp*sin (M2F*t))

# Use the "GfsOcean" model
1 0 GfsOcean GfsBox GfsGEdge {} {
    # Set the timestep to sthg small compared to the tidal period
    Time { dtmax = 100. }

    # Set the box size to 500 km
    PhysicalParams { L = 500e3 }

    # Use cartographic projection module
    GModule map

    # Use a Lambert conformal conic projection centered on 174 degrees
    # longitude and -40.8 degrees latitude. Rotate the domain 40
    # degrees counter-clockwise.
    MapProjection { lon = 174 lat = -40.8 angle = 40 }

    # Refine to six levels
    Refine 6

    # We want more accuracy in the projection than the default 1e-3
    ApproxProjectionParams { tolerance = 1e-6 nitermax = 10 }

    # Initialise tidal amplitudes
    Init {} {
        A_amp = AM2.gts
        B_amp = BM2.gts
    }

    # Bathymetry
    Solid bath.gts

    # Refine the coastline to 10 levels
    RefineSurface 10 bath.gts { twod = 1 }

    # Acceleration of gravity
    PhysicalParams { g = 9.81 }

    # Add Coriolis source term
    SourceCoriolis -1e-4

    # Bottom friction parameterisation:
    # Quadratic drag with friction coefficient of 4e-4.
    Init { istep = 1 } {
        U = U/(1. + dt*Velocity*4e-4/H)
        V = V/(1. + dt*Velocity*4e-4/H)
    }

    # Weak exponential filtering of the velocity field
    #    EventFilter { istep = 1 } U 40000
    #    EventFilter { istep = 1 } V 40000
    
    # After t=100000, starts on-the-fly harmonic decomposition of the pressure field...
    EventHarmonic { start = 100000 istep = 10 } P A B Z EP M2F
    # ... and of the velocity field
    EventHarmonic { start = 100000 istep = 10 } U AU BU ZU EU M2F
    EventHarmonic { start = 100000 istep = 10 } V AV BV ZV EV M2F

    # After t=100000, stops the simulation if the variations of the A0
    # harmonic component are less than 0.025 in 100 timesteps.
    EventStop { start = 100000 istep = 100 } A0 0.025

    OutputTime { istep = 1 } stderr
    OutputProjectionStats { istep = 1 } stderr

    # Output solution on standard output every 20 timesteps
    # for on-the-fly visualisation with GfsView.
    # Do not include the GTS file for the embedded surface to save bandwidth.
    OutputSimulation { istep = 20 } stdout { solid = 0 }

    # Output solution in file 'end.gfs' at the end of the simulation
    OutputSimulation { start = end } end.gfs

    # Output curves using gnuplot
    EventScript { start = end } {
        cat <<EOF | gnuplot
        set term postscript eps lw 3 solid 20 colour
        set output 'pv.eps'
        set xlabel 'Time (days)'
        set ylabel 'Elevation (metres) or Velocity (metres/s)'
        plot 'p' u (\$3/86400.):(\$9/9.81) w l t "Elevation", \
             'u' u (\$3/86400.):9 w l t "Velocity"
        set output 'a0.eps'
        set ylabel 'Maximum harmonic elevation (metres)'
        plot [1:]'a0' u (\$3/86400.):(\$9/9.81) w l t ""
EOF
    }

    OutputScalarNorm { istep = 1 } p { v = P }
    OutputScalarNorm { istep = 1 } u { v = Velocity }
    OutputScalarNorm { istep = 10 } a0 { v = sqrt(A0*A0 + B0*B0) }
}
GfsBox {
    # Use Flather boundary conditions on all boundaries
    left = Boundary {
        BcFlather U 0 H P M2(t)
    }
    right = Boundary {
        BcFlather U 0 H P M2(t)
    }
    top = Boundary {
        BcFlather V 0 H P M2(t)
    }
    bottom = Boundary {
        BcFlather V 0 H P M2(t)
    }

    # This is required for consistent free-surface fluxes
    front = Boundary
}