Creates a Parameter struct with following options and default values

T::DataType=Float32                 # number format

Tprog::DataType=T                   # number format for prognostic variables
Tcomm::DataType=Tprog               # number format for ghost-point copies

# DOMAIN RESOLUTION AND RATIO
nx::Int=100                         # number of grid cells in x-direction
Lx::Real=2000e3                     # length of the domain in x-direction [m]
L_ratio::Real=2                     # Domain aspect ratio of Lx/Ly

# PHYSICAL CONSTANTS
g::Real=10.                         # gravitational acceleration [m/s]
H::Real=500.                        # layer thickness at rest [m]
ρ::Real=1e3                         # water density [kg/m^3]
ϕ::Real=45.                         # central latitude of the domain (for coriolis) [°]
ω::Real=2π/(24*3600)                # Earth's angular frequency [s^-1]
R::Real=6.371e6                     # Earth's radius [m]

# WIND FORCING OPTIONS
wind_forcing_x::String="channel"    # "channel", "double_gyre", "shear","constant" or "none"
wind_forcing_y::String="constant"   # "channel", "double_gyre", "shear","constant" or "none"
Fx0::Real=0.12                      # wind stress strength [Pa] in x-direction
Fy0::Real=0.0                       # wind stress strength [Pa] in y-direction
seasonal_wind_x::Bool=false         # Change the wind stress with a sine of frequency ωFx,ωFy
seasonal_wind_y::Bool=false         # same for y-component
ωFx::Real=1.0                       # frequency [1/year] for x component
ωFy::Real=1.0                       # frequency [1/year] for y component

# BOTTOM TOPOGRAPHY OPTIONS
topography::String="ridge"          # "ridge", "seamount", "flat", "ridges", "bathtub"
topo_height::Real=50.               # height of seamount [m]
topo_width::Real=300e3              # horizontal scale [m] of the seamount

# SURFACE RELAXATION
surface_relax::Bool=false           # yes?
t_relax::Real=100.                  # time scale of the relaxation [days]
η_refh::Real=5.                     # height difference [m] of the interface relaxation profile
η_refw::Real=50e3                   # width [m] of the tangent used for the interface relaxation

# SURFACE FORCING
surface_forcing::Bool=false         # yes?
ωFη::Real=1.0                       # (annual) frequency [1/year]
A::Real=3e-5                        # Amplitude [m/s]
ϕk::Real=ϕ                          # Central latitude of Kelvin wave pumping
wk::Real=10e3                       # width [m] in y of Gaussian used for surface forcing

# TIME STEPPING OPTIONS
RKo::Int=4                          # Order of the RK time stepping scheme (3 or 4)
cfl::Real=1.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
Ndays::Real=10.0                    # number of days to integrate for
nstep_diff::Int=1                   # diffusive part every nstep_diff time steps.
nstep_advcor::Int=0                 # advection and coriolis update every nstep_advcor time steps.
                                    # 0 means it is included in every RK4 substep

# BOUNDARY CONDITION OPTIONS
bc::String="periodic"               # "periodic" or anything else for nonperiodic
α::Real=2.                          # lateral boundary condition parameter
                                    # 0 free-slip, 0<α<2 partial-slip, 2 no-slip

# MOMENTUM ADVECTION OPTIONS
adv_scheme::String="ArakawaHsu"     # "Sadourny" or "ArakawaHsu"
dynamics::String="nonlinear"        # "linear" or "nonlinear"

# BOTTOM FRICTION OPTIONS
bottom_drag::String="quadratic"     # "linear", "quadratic" or "none"
cD::Real=1e-5                       # bottom drag coefficient [dimensionless] for quadratic
τD::Real=300.                       # bottom drag coefficient [days] for linear

# DIFFUSION OPTIONS
diffusion::String="constant"        # "Smagorinsky" or "constant", biharmonic in both cases
νB::Real=500.0                      # [m^2/s] scaling constant for constant biharmonic diffusion
cSmag::Real=0.15                    # Smagorinsky coefficient [dimensionless]

# TRACER ADVECTION
tracer_advection::Bool=true         # yes?
tracer_relaxation::Bool=false       # yes?
tracer_consumption::Bool=false      # yes?
tracer_pumping::Bool=false          # yes?
injection_region::String="west"     # "west", "south", "rect" or "flat"
sst_initial::String="south"         # "west", "south", "rect", "flat" or "restart"
sst_rect_coords::Array{Float64,1}=[0.,0.15,0.,1.0]
                                    # (x0,x1,y0,y1) are the size of the rectangle in [0,1]
Uadv::Real=0.25                     # Velocity scale [m/s] for tracer advection
SSTmax::Real=1.                     # tracer (sea surface temperature) max for restoring
SSTmin::Real=0.                     # tracer (sea surface temperature) min for restoring
τSST::Real=500.                     # tracer restoring time scale [days]
jSST::Real=365.                     # tracer consumption [days]
SST_λ0::Real=222e3                  # [m] transition position of relaxation timescale
SST_λs::Real=111e3                  # [m] transition width of relaxation timescale
SST_γ0::Real=8.35                   # [days] injection time scale
SSTw::Real=5e3                      # width [m] of the tangent used for the IC and interface relaxation
SSTϕ::Real=0.5                      # latitude/longitude fraction ∈ [0,1] of sst edge

# OUTPUT OPTIONS
output::Bool=false                  # netcdf output?
output_vars::Array{String,1}=["u","v","η","sst","q","ζ"]  # which variables to output? "du","dv","dη" also allowed
output_dt::Real=6                   # output time step [hours]
outpath::String=pwd()               # path to output folder

# INITIAL CONDITIONS
initial_cond::String="rest"         # "rest" or "ncfile" for restart from file
initpath::String=outpath            # folder where to pick the restart files from
init_run_id::Int=0                  # run id for restart from run number
init_starti::Int=-1                 # timestep to start from (-1 meaning last)
get_id_mode::String="continue"      # How to determine the run id: "continue" or "fill"
run_id::Int=-1                      # Output with a specific run id
init_interpolation::Bool=true       # Interpolate the initial conditions in case grids don't match?

u,v,η,sst = RunModel()

runs ShallowWaters with default parameters as defined in src/DefaultParameters.jl

Examples

julia> u,v,η,sst = RunModel(Float64,nx=200,output=true)

ArakawaHsu variables collected in a struct.

Generator function for ArakawaHsu VarCollection.

Bernoulli variables collected in a struct.

Generator function for Bernoulli VarCollection.

Bottomdrag variables collected in a struct.

Generator function for Bottomdrag VarCollection.

Generator function for the mutable struct Constants.

Generator function for the Grid struct.

Laplace variables collected in a struct.

Generator function for Laplace VarCollection.

Generator function for "empty" NcFiles struct.

Generator function for NcFiles struct, creating the underlying netCDF files.

P = ProgVars{T}(u,v,η,sst)

Struct containing the prognostic variables u,v,η and sst.

Zero generator function for Grid G as argument.

Runge Kutta time stepping scheme diagnostic cariables collected in a struct.

Generator function for RungeKutta VarCollection.

SemiLagrange variables collected in a struct.

Generator function for SemiLagrange VarCollection.

Smagorinsky variables collected in a struct.

Generator function for Smagorinsky VarCollection.

Tendencies collected in a struct.

Generator function for Tendencies VarCollection.

VolumeFluxes collected in a struct.

Generator function for VolumeFluxes VarCollection.

Vorticity variables collected in a struct.

Generator function for Vorticity VarCollection.

Convert function for two arrays, X1, X2, in case their eltypes differ. Convert every element from X1 and store it in X2.

Convert function for two arrays, X1, X2, in case their eltypes are identical. Just pass X1, such that X2 is pointed to the same place in memory.

Linear combination α of potential voriticity q according to the energy and enstrophy conserving scheme of Arakawa and Hsu, 1990

Linear combination δ of potential voriticity q according to the energy and enstrophy conserving scheme of Arakawa and Hsu, 1990

Linear combination γ of potential voriticity q according to the energy and enstrophy conserving scheme of Arakawa and Hsu, 1990

Linear combination β of potential voriticity q according to the energy and enstrophy conserving scheme of Arakawa and Hsu, 1990

Returns the constant forcing matrices Fx,Fy that vary only meriodionally/zonally as a cosine with strongest forcing in the middle and vanishing forcing at boundaries.

Returns constant in in space forcing matrices Fx,Fy.

Returns the constant forcing matrices Fx,Fy that vary only meriodionally/zonally with a superposition of sin & cos for a double gyre circulation. See Cooper&Zanna 2015 or Kloewer et al 2018.

Returns a matrix of constant water depth H.

Time evolution of forcing.

Returns a reference state for Newtonian cooling/surface relaxation shaped as a hyperbolic tangent to force the continuity equation.

Linear interpolation of a variable u in the x-direction. m,n = size(ux) must be m+1,n = size(u).

Bilinear interpolation a variable u in x and y-direction. m,n = size(uxy) must be m+1,n+1 = size(u).

Linear interpolation a variable u in the y-direction. m,n = size(uy) must be m,n+1 = size(u).

Returns Kelvin wave pumping forcing of the continuity equation.

Returns constant in in space forcing matrices Fx,Fy.

Potential vorticity calculated as q = (f + ∂v/∂x - ∂u/∂y)/h.

Advection of potential vorticity qhv,qhu as in Arakawa and Hsu, 1990 Energy and enstrophy conserving (in the limit of non-divergent mass flux) scheme with τ = 0.

Advection of potential vorticity qhv,qhu as in Sadourny, 1975 enstrophy conserving scheme.

Transit function to call the specified advection scheme.

Returns a matrix of water depth for the whole domain that contains a meridional Gaussian ridge in the middle. Water depth, heigth and width of the ridge are adjusted with the constants waterdepth, topofeatheight and topofeat_width.

Same as Ridge() but for 3 ridges at 1/4,1/2,3/4 of the domain.

Returns a matrix of water depth for the whole domain that contains a Gaussian seamount in the middle. Water depth, heigth and width of the seamount are adjusted with the constants H, topofeatheight and topofeatwidth.

Returns the constant forcing matrices Fx,Fy that vary only meriodionally/zonally as a hyperbolic tangent with strongest shear in the middle.

Zonal mass flux U = uh.

Meridional mass flux V = vh.

Update u with bottom friction tendency (Bu,Bv) and biharmonic viscosity.

Extends the matrices u,v,η,sst with a halo of ghost points for boundary conditions.

Advection of sst/tracer based on the departure points xd,yd via bilinear interpolation. Departure points are clipped/wrapped to remain within the domain. Boundary conditions either periodic (wrap around behaviour) or no-flux (no gradient via clipping). Once the respective 4 surrounding grid points are found do bilinear interpolation on the unit square.

Update advective and Coriolis tendencies.

Add to a x multiplied with b. a += x*b

Solves the trajectory equation for a given time step dt and the velocity uv (this can be u or v). One function for three cases

(i) u is interpolated from u-grid with halo onto T-grid (ii) v is interpolated from v-grid with halo onto T-grid (iii) u or v already on the T-grid: All matrices have same size.

Uses relative grid nodes in the departure points rd, such that actually solved is rd = 0 - dt*uv. The indices i,j of rd determine the arrival grid point.

Bernoulli potential p = 1/2*(u² + v²) + gη.

Bilinear interpolation on (x,y) in the unit square [0,1]x[0,1]. The values at the corners are f00 = f(0,0), f01 = f(0,1), etc.

Transit function to call the specified bottom drag function.

Linear bottom drag computed as rB*(u,v). rB is negative and contains the grid spacing Δ as gradient operators are dimensionless.

Quadratic bottom drag Bu,Bv = cD/h * | u⃗ | * u⃗

c equals add a to x multiplied with b. c = a + x*b

Clips the relative lower-left corner index xyi (for both x or y indices) to remain within the domain. ij is the index of the underlying matrix. xy is the actual coordinate, mn (m or n) the size of the domain, and h1,h2 are the halo sizes (left/south and right/north).

Transit function to call the specified continuity function.

Continuity equation's right-hand side with time&space dependent forcing.

Continuity equation's right-hand side -∂x(uh) - ∂y(vh) without forcing.

Continuity equation's right-hand side with surface relaxation -∂x(uh) - ∂y(vh) + γ*(η_ref - η).

Coriolis parameter f [1/s] at latitude ϕ [°] given Earth's rotation ω [1/s].

Computes the departure point for semi-Lagrangian advection following Diamantakis, 2014. u,v are assumed to be the time averaged velocities over the previous advection time step. (Presumably need to be changed to 2nd order extrapolation in case the tracer is not passive)

Uses fixed-point iteration once to find the departure point.

Transit function to call the specified diffusion scheme.

Biharmonic diffusion operator with constant viscosity coefficient Viscosity = ν∇⁴ ⃗u. Although constant, the coefficient is actually inside Viscosity = ∇⋅ν∇∇² ⃗u.

Smagorinsky-like biharmonic viscosity Viscosity = ∇ ⋅ (cSmag Δ⁴ |D| ∇∇² ⃗u) The Δ⁴-scaling is omitted as gradient operators are dimensionless.

Estimates the total time the model integration will take.

Feedback function that calls duration estimate, nan_detection and progress.

Finalises the progress txt file.

Initialises the progress txt file.

Coriolis term f*u.

Coriolis term f*v.

Finds the first gap in a list of integers.

Checks output folders to determine a 4-digit run id number.

Decide on boundary condition P.bc which ghost point function to execute.

Decide on boundary condition P.bc which ghost point function to execute.

Copy ghost points for η from inside to the halo in the nonperiodic case.

Copy ghost points for η from inside to the halo in the periodic case.

Copy ghost points for u from inside to the halo in the nonperiodic case.

Copy ghost points for u from inside to the halo in the periodic case.

Copy ghost points for v from inside to the halo in the nonperiodic case.

Copy ghost points for v from inside to the halo in the periodic case.

Copy ghost points for η from inside to the halo in the nonperiodic case.

Copy ghost points for η from inside to the halo in the periodic case.

Interpolates the matrix uv into the matrix uvi, where xx, yy specify the coordinates as indices (including fraction. Interpolation only onto the inner entries of uvi. (They will be copied back later via the ghostpoint function). Two cases (i) u velocities: from u-grid with halo to T-grid (ii) v velocities: from v-grid with halo to T-grid.

Meter per 1 degree of latitude (or longitude at the equator).

Similar to the numpy meshgrid function: repeats x length(y)-times and vice versa. Returns two matrices xx,yy of same shape so that each row of xx is x and each column of yy is y.

Sum up the tendencies of the non-diffusive right-hand side for the u-component.

Sum up the tendencies of the non-diffusive right-hand side for the v-component.

Returns a boolean whether the prognostic variables contains a NaN.

Creates a netCDF file based on grid vectors x,y the variable name, its path, unit and long_name.

No bottom drag.

Closes netCDF and progress.txt files.

Preallocate the diagnostic variables and return them as matrices in structs.

Converts time step into percent for feedback.

Returns a human readable string representing seconds in terms of days, hours, minutes or seconds.

Cut off the halo from the prognostic variables.

Transit function to call either the rhslinear or the rhsnonlinear.

Tendencies du,dv,dη of

    ∂u/∂t = gv - g∂η/∂x + Fx
    ∂v/∂t = -fu - g∂η/∂y
    ∂η/∂t =  -∂(uH)/∂x - ∂(vH)/∂y + γ(η_ref-η) + Fηt*Fη,

the linear shallow water equations.

Tendencies du,dv,dη of

    ∂u/∂t = qhv - ∂(1/2*(u²+v²) + gη)/∂x + Fx
    ∂v/∂t = -qhu - ∂(1/2*(u²+v²) + gη)/∂y + Fy
    ∂η/∂t =  -∂(uh)/∂x - ∂(vh)/∂y + γ(η_ref-η) + Fηt*Fη

the nonlinear shallow water equations.

νSmag = cSmag * |D|, where deformation rate |D| = √((∂u/∂x - ∂v/∂y)^2 + (∂u/∂y + ∂v/∂x)^2). The grid spacing Δ is omitted here as the operators are dimensionless.

Squared velocities u²,v².

Biharmonic stress tensor ∇∇²(u,v) = (∂/∂x(∇²u), ∂/∂y(∇²u); ∂/∂x(∇²v), ∂/∂y(∇²v))

Layer thickness h obtained by adding sea surface height η to bottom height H.

Integrates ShallowWaters forward in time.

Tracer consumption via relaxation back to .

Biharmonic stress tensor times constant viscosity coefficient νB * ∇∇² ⃗u = (S11, S12; S21, S22)

Biharmonic stress tensor times Smagorinsky coefficient νSmag * ∇∇² ⃗u = (S11, S12; S21, S22).

Clips the relative lower-left corner index xyi (for both x or y indices) to remain within the domain. ij is the index of the underlying matrix. xy is the actual coordinate, mn (m or n) the size of the domain, and h is the halo size.

Helper function to create yy_q based on the boundary condition bc.

Coriolis parameter's derivative β wrt latitude [(ms)^-1] at latitude ϕ, given Earth's rotation ω [1/s] and radius R [m].

Calculates the 2nd order centred gradient in x-direction on any grid (u,v,T or q). The size of dudx must be m-1,n compared to m,n = size(u)

Calculates the 2nd order centred gradient in y-direction on any grid (u,v,T or q). The size of dudy must be m,n-1 compared to m,n = size(u).

∇² is the 2nd order centred Laplace-operator ∂/∂x^2 + ∂/∂y^2. The 1/Δ²-factor is omitted and moved into the viscosity coefficient.