GfsRiver

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 04:02, 28 August 2009
Popinet (Talk | contribs)
(Added description)
← Previous diff
Revision as of 20:15, 14 January 2010
Popinet (Talk | contribs)
(Added doc for optional parameters)
Next diff →
Line 10: Line 10:
The value of the acceleration of gravity is set using [[GfsPhysicalParams]] (default to 1). The value of the acceleration of gravity is set using [[GfsPhysicalParams]] (default to 1).
 +
 +The syntax in parameter files is
 +
 + [ [[GfsSimulation]] ] {
 + time_order = 2
 + dry = 1e-6
 + }
 + GfsBox {}
 + ...
 +
 +where the parameter block is optional. The optional parameters are:
 +
 +; <code>time_order</code>: the order of the temporal discretisation, either 1 or 2. The default is 2.
 +; <code>dry</code>: the depth threshold below which a cell is considered to be "dry". The default is 1e-6 times the box length scale (<code>L</code> parameter in [[GfsPhysicalParams]]).
<examples/> <examples/>

Revision as of 20:15, 14 January 2010

GfsRiver solves the Saint-Venant equations (also known as the non-linear shallow-water equations) in conservative form. Wetting and drying, hydraulic jumps etc... are handled appropriately by the solver.

The default variables are:

U
horizontal component of the flux (dimension L^2/T)
V
vertical component of the flux
P
fluid depth
Zb
bed elevation
H
Zb + P i.e. the elevation of the free surface

The value of the acceleration of gravity is set using GfsPhysicalParams (default to 1).

The syntax in parameter files is

[ GfsSimulation ] {
  time_order = 2
  dry = 1e-6
}
GfsBox {}
...

where the parameter block is optional. The optional parameters are:

time_order
the order of the temporal discretisation, either 1 or 2. The default is 2.
dry
the depth threshold below which a cell is considered to be "dry". The default is 1e-6 times the box length scale (L parameter in GfsPhysicalParams).

Examples

  • Dam break on complex topography
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
        PhysicalParams { L = 8. }
        
        # Define a 1D domain using cell masking
        RefineSurface 9 (y - 8.*(0.5 - 1./512.))
        InitMask {} (y < 8.*(0.5 - 1./512.))
    
        # Set the topography Zb and the initial water surface elevation P
        Init {} {
    	Zb = x*x/8.+cos(M_PI*x)/2.
    	P = {
    	    double p = x > 0. ? 0.35 : x < -2. ? 1.9 : -1.;
    	    return MAX (0., p - Zb);
    	}
        }
        PhysicalParams { g = 1. }
    
        # Use a first-order scheme rather than the default second-order
        # minmod limiter. This is just to add some numerical damping.
        AdvectionParams {
           # gradient = gfs_center_minmod_gradient
    	gradient = none
        }
    
        Time { end = 40 }
        OutputProgress { istep = 10 } stderr
        OutputScalarSum { istep = 10 } ke { v = (P > 0. ? U*U/P : 0.) }
        OutputScalarSum { istep = 10 } vol { v = P }
        OutputScalarNorm { istep = 10 } u { v = (P > 0. ? U/P : 0.) }
    
        # use gfsplot/gnuplot for online visualisation and generation of figures
        OutputSimulation { step = 0.3 } {
    	gfsplot "
              load 'dam.plot'
              set title sprintf('t = %4.1f', (t))
              plot [-4.:4.]'-' u (x):(Zb):(H) w filledcu lc 3, \
                           '-' u (x):(Zb) w l lw 4 lc 1 lt 1
              set term pngcairo size 800,600
              set output sprintf('sim-%04.1f.png', (t))
              replot
              set term wxt noraise
            "
        } { format = text }
    
        # Combine all the gif images into a gif animation using gifsicle
        EventScript { start = end } {
    	echo "\rcreating animation...            " > /dev/stderr
            sleep 10 # give a chance to gnuplot to catch up
    	for f in sim-*.png; do
    	    convert $f -trim +repage -bordercolor white \
    		-border 10 -resize 640x282! `basename $f .png`.gif && rm -f $f
    	done
    	gifsicle --colors 256 --optimize --delay 12 --loopcount=0 sim-*.gif > dam.gif && \
    	    rm -f sim-*.gif
        }
    }
    

  • Small amplitude solitary wave interacting with a parabolic hump
  • 2 1 GfsRiver GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
        Refine 8
        Init {} {
    	# Parabolic hump
    	Zb = 0.8*exp(-5.*(x - 0.9)*(x - 0.9) - 50.*(y - 0.5)*(y - 0.5))
    	# Initial free surface and perturbation
    	P = (0.05 < x && x < 0.15 ? 1.01 : 1) - Zb
        }
        PhysicalParams { g = 1 }
        AdvectionParams { cfl = 0.5 }
        AdaptGradient { istep = 1 } { 
    	cmax = 1e-4
    	cfactor = 2
    	maxlevel = 8
    	minlevel = 6
        } (P + Zb)
        Time { end = 1.8 }
        OutputTime { istep = 10 } stderr
        OutputSimulation { istep = 10 } stdout
        EventScript { istep = 10 } { echo "Save stdout { width = 640 height = 480 }" }
        OutputSimulation { start = 0.6 step = 0.3 } sim-%g.gfs
        EventScript { start = end } {
    	for i in 0.6 0.9 1.2 1.5 1.8; do
    	    echo "Save stdout { format = EPS line_width = 0.2 }" | \
    		gfsview-batch2D sim-$i.gfs isolines.gfv > iso-$i.eps
    	    echo "Save stdout { format = EPS line_width = 0.2 }" | \
    		gfsview-batch2D sim-$i.gfs cells.gfv > cells-$i.eps
    	done
    	echo "Save stdout { width = 1280 height = 960 }" | \
    	    gfsview-batch2D sim-0.9.gfs hump.gfv | convert ppm:- hump.eps
        }
    }
    

  • Shock reflection by a circular cylinder
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
        Time { end = 0.3 }
    
        PhysicalParams { L = 5 g = 9.81 }
    
        # We use a sphere knowing that in 2D the resulting object will be
        # a cross-section of the sphere at z = 0 i.e. a cylinder of radius
        # 0.5
        Solid (sphere(-0.5,0.,0.,0.5))
    
        Init {} {
    	# Initial shock
    	P = (x <= -1. ? 3.505271526 : 1.)
    	U = (x <= -1. ? 22.049341608 : 0.)
        }
    
        Refine 9
        AdaptGradient { istep = 1 } {
            cmax = 0.1
            cfactor = 2
            maxlevel = 9
        } P
    
        OutputTime { istep = 10 } stderr
        OutputSimulation { istep = 10 } stdout
    
        GModule gfsview
        OutputView { step = 0.0025 } {
    	ppm2mpeg -b 3600K > depth.mpg
        } { width = 400 height = 400 } depth.gfv
        OutputView { start = end } { 
    	convert ppm:- depth.eps 
        } { width = 400 height = 400 } depth.gfv
        OutputView { start = end } { 
    	convert ppm:- mesh.eps 
        } { width = 400 height = 400 } mesh.gfv
    }
    

  • Tsunami runup onto a complex three-dimensional beach
  • 2 1 GfsRiver GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
        # the domain is 3.402 m X 6.804 m
        # units for time are seconds
        PhysicalParams { L = 3.402 g = 9.81 }
        Refine 6
    
        # maintain the Zb1 variable using the GTS surface
        VariableFunction Zb1 bathy.gts
    
        # the initial water level is at z = 0, so the depth P is...
        Init {} {
            P = MAX (0., -Zb1)
        }
    
        # use a Sweby limiter rather than the default minmod which is too
        # dissipative
        AdvectionParams { gradient = gfs_center_sweby_gradient }
    
        # adapt down to 9 levels based on the slope of the (wet)
        # free-surface and with a tolerance of 1 mm
        AdaptGradient { istep = 1 } {
            cmax = 1e-3
            cfactor = 2
            maxlevel = 9
            minlevel = 6
        } (P < DRY ? 0. : P + Zb)
    
        # at each timestep
        Init { istep = 1 } {
    	# Add a "shelf" to simulate the wall on the right-hand-side boundary
    	Zb = (x > 5.448 ? 0.13535 : Zb1)
            # read in the experimental timeseries in variable 'input'
    	input = input.cgd
            # implicit quadratic bottom friction with coefficient 1e-3
    	U = (P > DRY ? U/(1. + dt*1e-3*Velocity/P) : 0.)
    	V = (P > DRY ? V/(1. + dt*1e-3*Velocity/P) : 0.)
            P = (P > DRY ? P : 0.)
        }
    
        Time { end = 22.5 }
        OutputTime { istep = 10 } stderr
        OutputTiming { start = end } stderr
        OutputSimulation { step = 1 } stdout
        OutputSimulation { step = 1 } sim-%g.gfs
        EventScript { step = 1 } { gzip -f sim-*.gfs }
    
        # output data at probe locations
        OutputLocation { istep = 1 } input 1e-3 1.7 0
        OutputLocation { istep = 1 } p5 4.521 1.196 0
        OutputLocation { istep = 1 } p7 4.521 1.696 0
        OutputLocation { istep = 1 } p9 4.521 2.196 0
    
        # kinetic energy
        OutputScalarSum { istep = 1 } ke { v = (P > 0. ? Velocity2*P : 0.) }
        # free-surface elevation
        OutputScalarStats { istep = 1 } p { v = (Zb > 0. ? P : P + Zb) }
        OutputScalarNorm { istep = 1 } u { v = Velocity }
        OutputTime { istep = 1 } balance
        OutputBalance { istep = 1 } balance
    
        # generate movies
        GModule gfsview
        OutputView { start = 9 step = 0.0416 } { ppm2mpeg -s 640x480 > monai.mpg } { 
    	width = 1280 height = 960
        } 3D.gfv
        OutputView { start = 14.63 end = 19.5 step = 0.033333333 } { 
    	ppm2mpeg -s 400x600 > overhead.mpg 
        } { 
    	width = 800 height = 1200
        } overhead.gfv
    } {
        dry = DRY
    }
    

  • The 2004 Indian Ocean tsunami
  • 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
    }
    

  • Poiseuille flow with multilayer Saint-Venant
  • 1 1 GfsRiver GfsBox GfsGEdge {} {
        Layers NL
        PhysicalParams { L = 1. }
        Init {} { P = 0.5 }
        Source U A*P/NL
        EventStop { istep = 1 } U 1e-8 DU
    #    OutputTime { istep = 1 } stderr
        Time { iend = 20 }
        OutputSimulation { start = end } {
    	awk '{
      	  nl = NL;
              if ($1 == "#") {
                for (i = 2; i <= NF; i++) {
                  split($i,a,":")
                  if (a[2] == "U0")
                    start = a[1];
                }
              }
              else {
                dz = $4/nl;
                emax = 0.;
      	    for (i = 0; i < nl; i++) {
                  z = dz*(0.5+i)
                  u = A/(2.*NU)*(1./4 - (0.5 - z)*(0.5 - z))
                  eu = u - $(start+i)/dz
                  if (eu < 0.)
                    eu = -eu;
                  if (eu > emax)
                    emax = eu;
                }
                print nl,emax
              }
            }'
        } { format = text }
    } {
        # vertical viscosity
        nu = NU
    }
    

  • Multi-layer Saint-Venant solver
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
        Layers NL
        PhysicalParams { L = RATIO g = 100./RE }
        Refine LEVEL
        InitMask {} (y < RATIO*(0.5 - 1./pow(2,LEVEL)))
    
        Init {} { P = 1 }
        Time { end = 1000 }
        EventStop { step = 1 } U0 1e-9
    #    OutputTime { step = 1 } stderr
    #    OutputSimulation { start = end } end-RATIO-RE.txt { format = text }
    #    OutputSimulation { start = end } end-RATIO-RE.gfs
        OutputSimulation { start = end } {
    	awk '{
      	  nl = NL;
              if ($1 == "#") {
                for (i = 2; i <= NF; i++) {
                  split($i,a,":")
                  if (a[2] == "U0")
                    start = a[1];
                }
              }
              else if ($1 == RATIO/2**(LEVEL + 1)) {
                dz = $4/nl;
      	    for (i = 0; i < nl; i++) {
                  z = dz*(0.5+i)
                  print $1,$2,z,$(start+i)/dz
                }
              }
            }' > uprof-RATIO-RE-NL
        } { format = text }
    } {
        # vertical viscosity
        nu = 1./RE
        # Neumann condition at the surface
        dut = 1.
    }
    

  • Wind-driven stratified lake
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
        Layers NL
        Refine LEVEL
        InitMask {} (y < RATIO*(0.5 - 1./pow(2,LEVEL)))
        VariableTracer DRHO
        PhysicalParams {
    	L = RATIO g = 1./(ALPHA*RE)
    	alpha = 1./(1. + DRHO)
        }
        Init {} {
    	P = 1
    	DRHO = (z > 0.75 ? 0 : ALPHA/ALPHAT)*P/NL
        }
        Time { end = 100 }
        OutputSimulation { start = end } end-NL.txt { format = text }
        OutputSimulation { start = end } {
    	awk '{
      	  nl = NL;
              if ($1 == "#") {
                for (i = 2; i <= NF; i++) {
                  split($i,a,":")
                  if (a[2] == "U0")
                    start = a[1];
                }
              }
              else if ($1 == RATIO/2**(LEVEL + 1)) {
                dz = $4/nl;
      	    for (i = 0; i < nl; i++) {
                  z = dz*(0.5+i)
                  print $1,$2,z,$(start+i)/dz
                }
              }
            }' > uprof-NL
        } { format = text }
    } {
        # vertical viscosity
        nu = 1./RE
        # Neumann condition at the surface
        dut = 1.
    }
    

  • Geostrophic adjustment with Saint-Venant
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
      Time { iend = 1580 dtmax = 1000 }
      Refine 6
      Global {
          #define L0 1000e3
          #define H0 1000
          #define G 0.01
          #define R0 100e3
          #define ETA0 599.5
          #define F0 1.0285e-4
      }
      AdvectionParams { gradient = gfs_center_gradient }
      PhysicalParams { L = L0 g = G }
      Init {} {
        # e-folding radius = 100 km
        # umax = 0.5 m/s = sqrt(200)*exp(-1/2)
        P = (H0 + ETA0*exp (-(x*x + y*y)/(R0*R0)))
        U = 2.*G*ETA0*y/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))*P
        V = - 2.*G*ETA0*x/(F0*R0*R0)*exp (-(x*x + y*y)/(R0*R0))*P
      }
      SourceCoriolis F0
     
      OutputErrorNorm { istep = 1 } { 
          awk '{print $3/86400. " " $9; fflush (stdout); }' > e
      } { v = P } {
        s = (H0 + ETA0*exp (-(x*x + y*y)/(R0*R0)))
        v = E
      }
      GModule gfsview
      OutputView {  istart = 100 iend = 500 istep = 100 } error-%ld.eps { format = EPS } geo.gfv
      OutputView {  istart = 1500 } error-%ld.eps { format = EPS } geo.gfv
      EventScript { start = end } {
        cat <<EOF | gnuplot
        set term postscript eps lw 3 color solid 20
        set output 'geo_error.eps'
        set xlabel 'Time (days)'
        set ylabel 'Maximum error on surface height (m)'
        plot 'e.ref' t 'ref' w l, 'e' t '' w l
    EOF
      }
    }
    

  • Oscillations in a parabolic container
  • 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
    }
    

  • Parabolic container with embedded solid
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
    
        # Analytical solution, see Sampson, Easton, Singh, 2006
        Global {
    	static gdouble Psi (double x, double y, 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))*R(x,y);
    	}
        }
        
        PhysicalParams { L = 10000 }
        Refine LEVEL
        Solid ({
    	    double line1 = SLOPE*x - y + 28000./pow(2.,LEVEL);
    	    double line2 = - SLOPE*x + y + 28000./pow(2.,LEVEL);
    	    return union (line1, line2);
    	    })
        Init {} {
    	Zb = h0*pow(R(x,y), 2.)/(a*a)
    	P = MAX (0., h0 + Psi (x, y, 0.) - Zb)
        }
        Init { istep = 1 } {
    	Pt = MAX (0., h0 + Psi (x, y, t) - Zb)
        }
        PhysicalParams { g = G }
        SourceCoriolis 0 tau
        Time { end = 6000 }
        OutputSimulation { start = 1500 } { 
    	awk '
            function atan(x)
            {
              return atan2(x,1.);
            }
            function pow(x,y)
            {
              return x**y;
            }
            {
              if ($1 == "#")
                print $0;
              else {
                printf ("%g", R($1,$2));
                for (i = 2; i <= NF; i++)
                  printf (" %s", $i);
                printf ("\n");
              }
            }' > sim-LEVEL-1500.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*cos (ANGLE) + V*sin(ANGLE) }
        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
    }
    

  • Transcritical flow over a bump
  • 1 0 GfsRiver GfsBox GfsGEdge { x = 0.5 } {
        PhysicalParams { L = 25 g = 9.81 }
        Refine LEVEL
        InitMask {} (y < 25.*(0.5 - 1./pow(2,LEVEL)))
        Init {} {
    	Zb = MAX(0., 0.2 - 0.05*(x - 10.)*(x - 10.))
    	P = 0.33 - Zb
        }
        AdvectionParams { 
    	# Sweby seems to work better than minmod for this test case
    	gradient = gfs_center_sweby_gradient 
        }
        EventStop { step = 1 } H 1e-5
        Time { end = 1000 }
    #    OutputTime { istep = 1 } stderr
        OutputSimulation { start = end } {
    	gfsplot "
              set term postscript eps color solid lw 2 16
              set output 'h-LEVEL.eps'
              set xlabel 'x'
              set ylabel 'z'
              plot [5:15]'swashes' u 1:4 t 'topography' w l, \
                         'swashes' u 1:6 t 'free surface (analytical)' w l, \
                         '-' u (x):(H) w p t 'free surface (numerical)'
            "
        } { format = text }
        OutputErrorNorm { start = end } { awk '{ print LEVEL,$5,$7,$9 }' } { v = H } {
    	s = href.cgd
    	v = EH
        }
        OutputSimulation { start = end } end-LEVEL.txt { format = text }
    }
    

  • Transcritical flow with multiple layers
  • 1 0 GfsRiver GfsBox GfsGEdge { x = 0.5 } {
        Global {
    	#define LENGTH 21.
    	#define WIDTH 5.75
    	#define HEIGHT 0.2
    	#define Q 1.
    	#define HE 0.6
    	#define G 9.81
        }
        Layers NL
        PhysicalParams { L = LENGTH g = G }
        Refine LEVEL
        InitMask {} (y < LENGTH*(0.5 - 1./pow(2,LEVEL)))
        Init {} {
    	Zb = MAX(0., HEIGHT*(1. - 1./(WIDTH/2.)/(WIDTH/2.)*(x - 10.)*(x - 10.)))
    	P = 0.6 - Zb
        }
        AdvectionParams { 
    	# Sweby seems to oscillate in this case
    #	gradient = gfs_center_sweby_gradient 
        }
        # Stop when stationary solution is reached
        EventStop { step = 1 } H 1e-3
        Time { end = 1000 }
        # uncomment this for on-the-fly animation
        # OutputSimulation { step = 0.1 } {
        # 	gfsplot "
        #       set term wxt noraise
        #       set xlabel 'x'
        #       set ylabel 'z'
        #       set xrange [0:21]
        #       set xtics 0,2,20
        #       set key bottom right
        #       plot '-' u (x):(Zb) t 'topography' w l, \
        #            '-' u (x):(H) w p t 'free surface'
        #     "
        # } { format = text }
        OutputLocation { start = end } { awk -v nl=NL -f uprof.awk > uprof-10-LEVEL-NL } 10. 10.49 0
        OutputLocation { start = end } { awk -v nl=NL -f uprof.awk > uprof-15-LEVEL-NL } 15. 10.49 0
        OutputLocation { start = end } { awk -v nl=NL -f uprof.awk > uprof-20-LEVEL-NL } 20. 10.49 0
        OutputSimulation { start = end } {
    	awk '{ if ($1 != "#") print $1,$8; }' > prof-LEVEL-NL
        } { format = text }
        OutputSimulation { start = end } end-LEVEL-NL { format = text }
    } {
        # vertical viscosity
        nu = 0.01
        # Gauckler-Manning-Strickler bottom friction, G = 9.81 m/s^2,
        # Strickler coefficient S = 25 m^(1/3)/s.
        k = (G/(25.*25.*pow(P,1./3.)))*Velocity
    }
    

  • Tsunami runup onto a plane beach
  • 1 0 GfsRiver GfsBox GfsGEdge { x = 0.3333333 } {
        PhysicalParams { L = 60000 }
    
        # Set a solid boundary close to the top boundary to limit the
        # domain width to less than one cell (i.e. a 1D domain)
        RefineSolid (MINLEVEL + (LEVEL - MINLEVEL)*(1. - x/50000.))
        Solid (y - 60000.*(0.5 - 0.99/pow(2,LEVEL)))
    
        # Set the topography Zb and the initial water surface elevation P
        Init {} {
    	Zb = -x/10.
    	P = init.cgd
    	P = MAX (0., P - Zb)
        }
        PhysicalParams { g = 9.81 }
    
        AdvectionParams {
    	gradient = gfs_center_sweby_gradient
        }
    
        # Force the flow to stay 1D
        Init { istep = 10 } { V = 0 }
    
        Time { end = 220 }
    
        # OutputTime { istep = 10 } stderr
        OutputSimulation { start = 0 } sim-LEVEL-%g.txt { format = text }
        OutputSimulation { start = 160 } sim-LEVEL-%g.txt { format = text }
        OutputSimulation { start = 175 } sim-LEVEL-%g.txt { format = text }
        OutputSimulation { start = 220 } sim-LEVEL-%g.txt { format = text }
    } {
        dry = DRY
        time_order = 2
    }
    

  • Lake-at-rest balance in an inclined domain with cut cells
  • 1 0 GfsRiver GfsBox GfsGEdge { x = 0.5 } {
        PhysicalParams { L = 10 g = 9.81 }
    
        Refine 5
        RefineSolid 7
        Solid (sphere(5.,0.,0.,2.0))
    
        Init {} {
    	Zb = x/10.
    	P = 12.1/(100. - 4.*M_PI)
        }
    
        Time { end = 200 }
    #    OutputSimulation { istep = 10 } stdout
        SourceCoriolis 0 1.0e-01
        OutputScalarNorm { start = end } u { v = U }
        OutputErrorNorm { start = end } ep { v = P } {
    	s = MAX(0, 0.5 - x/10.)
    	unbiased = 1
    	relative = 1
        }
        GModule gfsview
        OutputView { start = end } still.eps { format = EPS } still.gfv
        EventScript { start = end } {
    	status=0
    	if awk '{if ($9 > 1e-5) { print "u: " $9 > "/dev/stderr"; exit (1); }}' < u ; then :
            else
                status=$GFS_STOP;
            fi
    	if awk '{if ($9 > 4e-3) { print "ep: " $9 > "/dev/stderr"; exit (1); }}' < ep ; then :
            else
                status=$GFS_STOP;
            fi
    	exit $status
        }
    }
    GfsBox {
        left = Boundary
        right = Boundary
        top = Boundary
        bottom = Boundary
    }
    

  • Lake-at-rest balance in an inclined domain with bipolar metric
  • 8 8 GfsRiver GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
        PhysicalParams { g = 9.81 }
    
        # Bipolar coordinates. see http://en.wikipedia.org/wiki/Bipolar_coordinates
        Metric M {
    	x = 10.*sinh(rx/2. + 1.)/(cosh(rx/2. + 1.) - cos (M_PI*ry/4.))
    	y = 10.*sin(M_PI*ry/4.)/(cosh(rx/2. + 1.) - cos (M_PI*ry/4.))
        }
        Refine 4
    
        Init {} {
    	Zb = x/10.
    	P = 12.1/(100. - 4.*M_PI)
        }
    
        Time { end = 200 }
    #    OutputSimulation { istep = 10 } stdout
        SourceCoriolis 0 1.0e-01
        OutputScalarNorm { start = end } u { v = U }
        OutputErrorNorm { start = end } ep { v = P } {
    	s = MAX(0, 1.26 - x/10.)
    	unbiased = 1
    	relative = 1
        }
        GModule gfsview
        OutputScalarSum { start = 0 } vol { v = Zb }
        OutputView { start = end } p.gnu { format = Gnuplot } p.gfv
        OutputView { start = end } mesh.gnu { format = Gnuplot } mesh.gfv
        EventScript { start = end } {
    	status=0
    	if gnuplot <<EOF; then :
            set term postscript eps lw 2 solid color
            set output 'still.eps'
            unset key
            unset logscale
            unset border
            unset xtics
            unset ytics
            unset xlabel
            unset ylabel
            set size ratio -1
            plot 'mesh.gnu' w l lc 0, 'p.gnu' w l lc 1
    EOF
    	else
    	    status=$GFS_STOP;
    	fi
    	if awk '{if ($9 > 8e-7) { print "u: " $9 > "/dev/stderr"; exit (1); }}' < u ; then :
            else
                status=$GFS_STOP;
            fi
    	if awk '{if ($9 > 3e-4) { print "ep: " $9 > "/dev/stderr"; exit (1); }}' < ep ; then :
            else
                status=$GFS_STOP;
            fi
    	exit $status
        }
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    GfsBox {
        left = Boundary
        right = Boundary
    }
    1 2 top
    2 3 top
    3 4 top
    4 5 top
    5 6 top
    6 7 top
    7 8 top
    8 1 top
    

  • Terrain reconstruction
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
        PhysicalParams { L = 8 }
        GModule terrain
        RefineTerrain LEVEL H {
    	basename = terrain,terrain-high
        } TRUE
        VariableTerrain T {
    	basename = terrain,terrain-high
        }
        Time { end = 0 }
        OutputErrorNorm { start = end } { awk '{ print LEVEL,$5,$7,$9 }' >> error-t } { v = T } {
    	s = {
    	    double r = sqrt (x*x + y*y);
    	    return r*r/8.+cos(M_PI*r)/2.;
    	}
    	w = (fabs (x) < 3.8 && fabs (y) < 3.8)
    	relative = 1
    	v = Et
        }
        OutputErrorNorm { start = end } { awk '{ print LEVEL,$5,$7,$9 }' >> error-h } { v = H0 } {
    	s = {
    	    double r = sqrt (x*x + y*y);
    	    return r*r/8.+cos(M_PI*r)/2.;
    	}
    	w = (fabs (x) < 3.8 && fabs (y) < 3.8)
    	relative = 1
    	v = Eh
        }
        OutputSimulation { start = end } {
    	awk '{
               if ($1 != "#") {
                 r = sqrt ($1*$1 + $2*$2);
                 print r,$4,$5,$6,$7;
               }
            }' | sort -n -k 1,2 > end-LEVEL.txt
        } { format = text variables = T,H0,Et,Eh }
        OutputSimulation { start = end } end-LEVEL.gfs
    }
    

  • Circular dam break on a sphere
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
        PhysicalParams { L = LENGTH }
        MetricLonLat M 1.
        Refine 8
        InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
        Init {} { P = 0.2 + 1.8*P/LENGTH }
        Time { end = 0.9 }
        OutputTime { istep = 10 } stderr
        GModule gfsview
        OutputView { start = 0.3 step = 0.3 } isolines-%g.eps { 
    	format = EPS line_width = 0.5 
        } isolines.gfv
        OutputSimulation { step = 0.3 } sim-%g.txt { variables = U,V,P format = text }
        OutputScalarSum { istep = 1 } sp { v = P }
    #    OutputSimulation { istep = 10 } stdout
        EventScript { start = end } {
    	status=0
    	for i in 0.3 0.6 0.9; do
    	    if awk 'BEGIN { n1 = 0; pi = 3.14159265359 }{
                  if ($1 != "#") {
                    c = cos($1*pi/180.)*cos($2*pi/180.);
                    d = atan2(sqrt(1. - c*c),c)*180./pi
                    i = int(d*2.)
                    x[i] += d
                    y[i] += $6
                    n[i]++
                    x1[n1] = d;
                    y1[n1++] = $6;
                  }
                }END {
                  for (i = 0; i <= 180; i++)
                    if (n[i] > 0)
                      print x[i]/n[i], y[i]/n[i], n[i];
                  sum = 0.;
                  for (i = 0; i < n1; i++) {
                    j = int(x1[i]*2.)
                    d = y1[i] - y[j]/n[j];
                    sum += d*d;
                  }
                  scatter = sqrt(sum/n1);
                  if (scatter > 0.013) {
                    print scatter > "/dev/stderr";
                    exit (1);
                  }
                }' < sim-$i.txt > prof-$i.txt ; then :
    	    else
    		status=$GFS_STOP;
    	    fi
    
    	    cat <<EOF | gnuplot
    rdist(x,y)=acos(cos(x*pi/180.)*cos(y*pi/180.))*180./pi
    cdist(x,y)=sqrt(x*x+y*y)
    set xlabel 'Angular distance (degree)'
    set ylabel 'Surface height'
    set xtics 0,22.5,90
    set ytics 0,0.25,0.75
    set term postscript eps color lw 2 solid 20
    set output 'p-$i.eps'
    plot [0:90][0:0.75]'sim-$i.txt' u (rdist(\$1,\$2)):6 ps 0.25 pt 6 t '', 'prof-$i.txt' w l lw 2 t ''
    EOF
    
    	done
    	exit $status
        }
    }
    

  • Circular dam break on a rotating sphere
  • 1 0 GfsRiver GfsBox GfsGEdge {} {
        PhysicalParams { L = LENGTH }
        MetricLonLat M 1.
        Refine 8
        InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
        Init {} { P = 0.2 + 1.8*P/LENGTH }
        Time { end = 1.4 }
        SourceCoriolis 10.*sin(y*M_PI/180.)
        OutputTime { istep = 10 } stderr
        GModule gfsview
        OutputView { start = 0.4 step = 0.4 } isolines-%g.eps { 
    	format = EPS line_width = 0.5 
        } isolines.gfv
        OutputSimulation { step = 0.4 } sim-%g.txt { variables = U,V,P format = text }
    #    OutputSimulation { istep = 10 } stdout
    }
    

  • Circular dam break on a ``cubed sphere''
  • 6 12 GfsRiver GfsBox GfsGEdge {} {
        PhysicalParams { L = 2.*M_PI/4. }
        MetricCubed M 8
        Refine 8
        InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
        Init {} { P = 0.2 + 1.8*P/LENGTH }
        Time { end = 1.5 }
        AdaptGradient { istep = 1 } { cmax = 1e-2 maxlevel = 8 } P
    #    OutputTime { istep = 10 } stderr
        GModule gfsview
        OutputView { start = 0.3 step = 0.3 } isolines-%g.gnu { 
    	format = Gnuplot
        } isolines.gfv
        OutputSimulation { step = 0.3 } sim-%g.txt { variables = U,V,P format = text }
        OutputScalarSum { istep = 1 } sp { v = P }
    #    OutputSimulation { istep = 10 } stdout
        EventScript { start = end } {
    	status=0
    	for i in 0.3 0.6 0.9 1.2 1.5; do
    	    if awk 'BEGIN { n1 = 0; pi = 3.14159265359 }{
                  if ($1 != "#") {
                    c = cos($1*pi/180.)*cos($2*pi/180.);
                    d = atan2(sqrt(1. - c*c),c)*180./pi
                    i = int(d*2.)
                    x[i] += d
                    y[i] += $6
                    n[i]++
                    x1[n1] = d;
                    y1[n1++] = $6;
                  }
                }END {
                  for (i = 0; i <= 180; i++)
                    if (n[i] > 0)
                      print x[i]/n[i], y[i]/n[i], n[i];
                  sum = 0.;
                  for (i = 0; i < n1; i++) {
                    j = int(x1[i]*2.)
                    d = y1[i] - y[j]/n[j];
                    sum += d*d;
                  }
                  scatter = sqrt(sum/n1);
                  if (scatter > 0.015) {
                    print scatter > "/dev/stderr";
                    exit (1);
                  }
                }' < sim-$i.txt > prof-$i.txt ; then :
    	    else
    		status=$GFS_STOP;
    	    fi
    
    	    gnuplot <<EOF
    set term postscript eps color lw 1 18
    set output 'isolines-$i.eps'
    set size ratio -1
    set xtics -90,30,90
    set ytics -90,30,90
    unset key
    plot [-90:90][-90:90]'isolines-$i.gnu' w l
    
    set term postscript eps color lw 2 solid 20
    set output 'p-$i.eps'
    rdist(x,y)=acos(cos(x*pi/180.)*cos(y*pi/180.))*180./pi
    cdist(x,y)=sqrt(x*x+y*y)
    set xlabel 'Angular distance (degree)'
    set ylabel 'Surface height'
    set xtics 0,22.5,90
    set ytics 0,0.25,0.75
    plot [0:90][0:0.75]'sim-$i.txt' u (rdist(\$1,\$2)):6 ps 0.25 pt 6 t '', 'prof-$i.txt' w l lw 2 t ''
    EOF
    
    	done
    	exit $status
        }
    }
    

  • Rossby--Haurwitz wave with Saint-Venant
  • 6 12 GfsRiver GfsBox GfsGEdge {} {
        Global {
    	#define AR 6.37122e6
    	#define N 4.
    	#define Umax 50.
    	#define M (Umax/(N*AR))
    	#define K M
    	
    	#define Omega 7.292e-5
    	#define DTR (M_PI/180.)
    	#define H0 8e3
    	#define G 9.81
    
    	// Williamson 1992, eq. (137)
    	static double u0 (double lambda, double phi) {
    	    double cosphi = cos (phi), sinphi = sin (phi);
    	    return AR*M*cosphi + AR*K*pow (cosphi, N - 1)*
    	    cos (N*lambda)*(N*sinphi*sinphi - cosphi*cosphi);
    	}
    
    	// Williamson 1992, eq. (138)
    	static double v0 (double lambda, double phi) {
    	    double cosphi = cos (phi), sinphi = sin (phi);
    	    return - AR*K*N*pow (cosphi, N - 1)*sinphi*sin (N*lambda);
    	}
    
            // Williamson 1992, eq. (139)
    	static double vorticity0 (double lambda, double phi) {
    	    return 2.*M*sin(phi) - K*sin(phi)*pow (cos (phi), N)*(N*N + 3.*N + 2.)*cos (N*lambda);
            }
    
            // Williamson 1992, eqs. (140-143)
    	static double p0 (double lambda, double phi, double t) {
    	    double nu = (N*(3. + N)*M - 2.*Omega)/((1. + N)*(2. + N));
    	    lambda -= nu*t;
    	    double cosphi = cos (phi);
    	    double Aphi = M/2.*(2.*Omega + M)*cosphi*cosphi + K*K/4.*pow (cosphi, 2.*N)*
    	      ((N + 1.)*cosphi*cosphi + 2.*N*N - N - 2. - 2.*N*N*pow(cosphi, -2.));
    	    double Bphi = 2.*(Omega + M)*K/((N + 1.)*(N + 2.))*pow(cosphi, N)*
                  (N*N + 2.*N + 2. - (N + 1.)*(N + 1.)*cosphi*cosphi);
    	    double Cphi = K*K/4.*pow(cosphi,2.*N)*((N + 1.)*cosphi*cosphi - (N + 2.));
    	    return AR*AR*(Aphi + Bphi*cos (N*lambda) + Cphi*cos (2.*N*lambda));
    	}
        }
        PhysicalParams { L = 2.*M_PI*AR/4. g = G }
        MetricCubed M LEVEL
        SourceCoriolis 2.*Omega*sin(y*DTR)
    
        Init {} { 
    	P = H0 + p0(x*DTR,y*DTR,t)/G
    	(U,V) = (u0(x*DTR,y*DTR)*P, v0(x*DTR,y*DTR)*P)
        }
    
        Refine LEVEL
    
        AdvectionParams { 
    	gradient = gfs_center_gradient 
    #	gradient = none
        }
    
        # ~24 days
        Time { end = 2073534 dtmax = 1e3 }
    #    OutputTime { istep = 10 } stderr
    
        OutputScalarNorm { istep = 10 } v-LEVEL { v = V }
        OutputScalarSum { istep = 10 } ec-LEVEL { v = Velocity2 }
    #    OutputScalarSum { istep = 10 } zeta-LEVEL { v = Vorticity }
        OutputScalarSum { istep = 10 } p-LEVEL { v = P }
        OutputErrorNorm { istep = 10 } eh-LEVEL { v = P } {
    	s = p0(x*DTR,y*DTR,t)/G
    	v = EH unbiased = 1 relative = 1
        }
        OutputSimulation { start = end } end-LEVEL.gfs
    #    OutputSimulation { istep = 10 } stdout
    #    OutputPPM { istep = 10 } eh.ppm { v = EH }
    #    OutputPPM { istep = 10 } p.ppm { v = P }
    
       GModule gfsview
       OutputView { start = end } ehp-LEVEL.gnu { format = Gnuplot } ehp.gfv
       OutputView { start = end } ehm-LEVEL.gnu { format = Gnuplot } ehm.gfv
       OutputView { start = end } h-LEVEL.gnu { format = Gnuplot } h.gfv
       OutputView { start = end } href-LEVEL.gnu { format = Gnuplot } href.gfv
    }
    

Personal tools
communication