# the radius of the Earth (in meters) 
Define RADIUS 6371220.
# the size of the domain
Define LENGTH 6000e3
# the domain is centered on 94E, 8N
Define LONGITUDE 94.
Define LATITUDE 8.
# anywhere with a fluid depth < 1 cm is considered dry
Define DRY 0.01
# 12 levels of refinement max i.e. a maximum resolution of 6000e3/2^12 = 1464 m
Define MAXLEVEL 12
# keep a coarse band, 0.04 wide, on all boundaries of the domain
# these act as "sponges" before the waves exit the domain
Define LEVEL (fabs (rx) < 0.46 && fabs (ry) < 0.46 ? MAXLEVEL : 5)

1 0 GfsRiver GfsBox GfsGEdge {} {
    PhysicalParams { 
	# length of the domain (m)
	L = LENGTH 
        # gravity is 9.81 m/s^2
	g = 9.81
	# from now on, units have been chosen to be metres and seconds
    }
    # use the longitude/latitude metric. x and y coordinates are now
    # longitude and latitude in degrees
    MetricLonLat M RADIUS
    # shift the origin to where we want it
    MapTransform { tx = LONGITUDE ty = LATITUDE }
    # 10 hours
    Time { end = 36000 }
    # use the terrain module for topography
    GModule terrain
    # use the okada module for initial deformations
    GModule okada
    Refine 5
    # maintain Zb as the terrain elevation defined by the samples in
    # ETOPO1
    # the 'etopo1' database needs to have been created beforehand and
    # be accessible within GFS_TERRAIN_PATH, see the 'terrain' module
    # documentation above
    VariableTerrain Zb {
        basename = etopo1
    } {
	# preserve lake-at-rest balance (but not volume) when
	# reconstructing bathymetry
	reconstruct = 1 
    }

    # the default minmod limiter is too dissipative for tsunamis
    AdvectionParams { gradient = gfs_center_sweby_gradient }

    # deformation from Grilli et al, 2007, Table 1
    # 5 segments triggered at different times

    # Segment 1
    Init {} { D = 0 }
    InitOkada D {
	x = 94.57 y = 3.83
	depth = 11.4857e3
	strike = 323 dip = 12 rake = 90
	length = 220e3 width = 130e3
	U = 18
    }
    # Initial water level is at z = D
    Init { start = 0 } { P = MAX (0., D - Zb) }

    # Segment 2
    EventList { start = 212 step = 6 end = 272 } {
	Init {} { D = 0 }
	InitOkada D {
	    x = 93.90 y = 5.22
	    depth = 11.4857e3
	    strike = 348 dip = 12 rake = 90
	    length = 150e3 width = 130e3
	    U = 23
	}
    }
    # make sure the deformation is well resolved
    AdaptGradient { start = 212 istep = 1 end = 272 } {
	cmax = 0.05 cfactor = 2
	minlevel = 5 maxlevel = LEVEL 
    } D
    # add it to the current elevation field (only if wet)
    Init { start = 272 } { P = (P < DRY ? P : MAX (0., P + D)) }

    # Segment 3
    EventList { start = 528 step = 6 end = 588 } {
	Init {} { D = 0 }
	InitOkada D {
	    x = 93.21 y = 7.41
	    depth = 12.525e3
	    strike = 338 dip = 12 rake = 90
	    length = 390e3 width = 120e3
	    U = 12
	}
    }
    # make sure the deformation is well resolved
    AdaptGradient { start = 528 istep = 1 end = 588 } {
	cmax = 0.05 cfactor = 2
	minlevel = 5 maxlevel = LEVEL 
    } D
    # add it to the current elevation field (only if wet)
    Init { start = 588 } { P = (P < DRY ? P : MAX (0., P + D)) }

    # Segment 4
    EventList { start = 853 step = 6 end = 913 } {
	Init {} { D = 0 }
	InitOkada D {
	    x = 92.60 y = 9.70
	    depth = 15.12419e3
	    strike = 356 dip = 12 rake = 90
	    length = 150e3 width = 95e3
	    U = 12
	}
    }
    # make sure the deformation is well resolved
    AdaptGradient { start = 853 istep = 1 end = 913 } {
	cmax = 0.05 cfactor = 2
	minlevel = 5 maxlevel = LEVEL 
    } D	
    # add it to the current elevation field (only if wet)
    Init { start = 913 } { P = (P < DRY ? P : MAX (0., P + D)) }

    # Segment 5
    EventList { start = 1213 step = 6 end = 1273 } {
	Init {} { D = 0 }
	InitOkada D {
	    x = 92.87 y = 11.70
	    depth = 15.12419e3
	    strike = 10 dip = 12 rake = 90
	    length = 350e3 width = 95e3
	    U = 12
	}
    }
    # make sure the deformation is well resolved
    AdaptGradient { start = 1213 istep = 1 end = 1273 } {
	cmax = 0.05 cfactor = 2
	minlevel = 5 maxlevel = LEVEL 
    } D	
    # add it to the current elevation field (only if wet)
    Init { start = 1273 } { P = (P < DRY ? P : MAX (0., P + D)) }

    # we don't want the arrival time to be interpolated from coarse
    # to fine, so we make it a "boolean" variable
    VariableBoolean Atime
    Init {} {
	P = MAX (0., D - Zb)
	Pmax = 0.
	Atime = -1.
    }

    Init { istep = 1 } {
	# elevation of the wet surface
	Hwet = (P > DRY ? H : NODATA)
	# maximum amplitude
	Pmax = (P > DRY && H > Pmax ? H : Pmax)
	# arrival time for an amplitude of +- 5 cm
	Atime = (Atime >= 0. ? Atime : P > DRY && fabs (H) > 5e-2 ? t : -1.)
	# implicit scheme for quadratic bottom friction with coefficient 1e-3
        U = P > DRY ? U/(1. + dt*Velocity*1e-3/P) : 0
        V = P > DRY ? V/(1. + dt*Velocity*1e-3/P) : 0
    }

    # adapt on gradient of wet surface elevation
    # with a tolerance of 5 cm
    AdaptGradient { istep = 1 } {
	cmax = 0.05 cfactor = 2
	minlevel = 5 maxlevel = LEVEL 
    } (P < DRY ? 0. : P + Zb)

    # also keep enough resolution to resolve the maximum elevation
    # field with a 5 cm discretisation error max
    AdaptError { istep = 1 } {
    	cmax = 0.05
    	minlevel = 5 maxlevel = LEVEL 
    } Pmax
    
    Global {
	#define RESOLUTION (LENGTH/pow (2, MAXLEVEL))
	#define close_enough(x,y) (distance(x,y,0.) < RESOLUTION)
    }
    # include a list of locations for which to output timeseries
    Include output.gfs

    # Jason profiles
    OutputLocation { start = 6900 } jason.txt jason.xy

    OutputTime { istep = 1 } stderr
    OutputBalance { istep = 1 } stderr
    OutputTiming { istep = 100 } stderr

    # save a snapshot every 100 iterations
    EventScript { istep = 100 } { rm -f snapshot-*.gfs }
    OutputSimulation { istep = 100 } snapshot-%ld.gfs

    # output solution on standard output every 20 timesteps
    # for on-the-fly visualisation with GfsView
    OutputSimulation { istep = 20 } stdout

    # ouput solution every 900 seconds (15 minutes)
    OutputSimulation { step = 900 } sim-%g.gfs
    EventScript { step = 900} { gzip -f sim-*.gfs }

    OutputScalarStats { istep = 1 } p { v = P }
    OutputScalarNorm { istep = 1 } u { v = Velocity }
    OutputScalarNorm { istep = 1 } U { v = U }
    OutputScalarNorm { istep = 1 } V { v = V }
    OutputScalarNorm { istep = 1 } hwet { v = Hwet }

    # animation of the wave elevation
    GModule gfsview
    OutputView { step = 60 } { ppm2mpeg > h.mpg } { 
	width = 800 height = 700
    } h.gfv
    OutputView { start = 7200 } { convert ppm:- eps2:h.eps } { 
	width = 800 height = 700
    } h.gfv
    # animation of the level of refinement
    OutputView { step = 60 } { ppm2mpeg > level.mpg } { 
	width = 800 height = 700
    } level.gfv
    OutputView { start = 7200 } { convert ppm:- eps2:level.eps } { 
	width = 800 height = 700 
    } level.gfv
    # animation of the velocity
    OutputPPM { step = 60 } { ppm2mpeg > velocity.mpg } {
	maxlevel = 10
	v = Velocity 
	min = 0 max = 0.25
    }
    # animation of the topography
    OutputPPM { step = 60 } { ppm2mpeg > zb.mpg } {
	maxlevel = 10
	v = Zb
    }
    # graphics
    OutputView { start = end } { convert ppm:- eps2:hmax.eps } { 
	width = 1600 height = 1600
    } hmax.gfv
    OutputView { start = end } { convert ppm:- eps2:hmax-detail.eps } { 
	width = 1600 height = 1600
    } hmax-detail.gfv
    EventScript { start = end } {
	# timeseries at specific locations
	gnuplot <<EOF
set term postscript eps color lw 2 size 5,2 18

set size ratio 0.4

set output 'hani.eps'
plot [3:5]'hanires.txt' u 1:(\$2/100.) w lp t 'observed' ps 2, \
           'hani.txt' u (\$3/3600.):7 w l t 'modelled'

set output 'male.eps'
plot [3:5]'maleres.txt' u 1:(\$2/100.) w lp t 'observed' ps 2, \
           'male.txt' u (\$3/3600.):7 w l t 'modelled'

set output 'gana.eps'
plot [3:5]'ganares.txt' u 1:(\$2/100.) w lp t 'observed' ps 2, \
           'gana.txt' u (\$3/3600.):7 w l t 'modelled'

set output 'dieg.eps'
plot [3:5]'diegres.txt' u 1:(\$2/100.) w lp t 'observed' ps 2, \
           'dieg.txt' u (\$3/3600.):7 w l t 'modelled'

set output 'colo.eps'
plot [2.5:4.5]'colores.txt' u 1:(\$2/100.) w lp t 'observed' ps 2, \
           'colo.txt' u (\$3/3600.):7 w l t 'modelled'

set output 'jason.eps'
set xlabel 'Latitude (deg)'
set ylabel 'Elevation (m)'
plot [-6:20][-1:1]'jasonres.txt' u 2:4 w l t 'observed', \
           'jason.txt' u 3:(\$3 > 19.2 ? 0 : \$8) w l t 'modelled'
EOF
    }
} {
    dry = DRY
}
GfsBox {
    # inflow/outflow on all boundaries with a reference water level at
    # z = 0
    top = Boundary {
	BcSubcritical V -Zb
    }
    bottom = Boundary {
	BcSubcritical V -Zb
    }
    right = Boundary {
	BcSubcritical U -Zb
    }
    left = Boundary {
	BcSubcritical U -Zb
    }
}