# Library

Documenting the public user interface.

## Oceananigans.jl

`Oceananigans.Oceananigans`

— ModuleMain module for `Oceananigans.jl`

– a Julia software for fast, friendly, flexible, data-driven, ocean-flavored fluid dynamics on CPUs and GPUs.

## Abstract operations

`Oceananigans.AbstractOperations.volume`

— Constant`volume = VolumeMetric()`

Instance of `VolumeMetric`

that generates `BinaryOperation`

s between `AbstractField`

s and their cell volumes. Summing this `BinaryOperation`

yields an integral of `AbstractField`

over the domain.

**Example**

```
julia> using Oceananigans
julia> using Oceananigans.AbstractOperations: volume
julia> c = CenterField(RectilinearGrid(size=(2, 2, 2), extent=(1, 2, 3)));
julia> c .= 1;
julia> c_dV = c * volume
BinaryOperation at (Center, Center, Center)
├── grid: 2×2×2 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
* at (Center, Center, Center)
├── 2×2×2 Field{Center, Center, Center} on RectilinearGrid on CPU
└── Vᶜᶜᶜ at (Center, Center, Center)
julia> c_dV[1, 1, 1]
0.75
julia> sum(c_dV)
6.0
```

`Oceananigans.AbstractOperations.Δz`

— Constant`Δz = ZSpacingMetric()`

Instance of `ZSpacingMetric`

that generates `BinaryOperation`

s between `AbstractField`

s and the vertical grid spacing evaluated at the same location as the `AbstractField`

.

`Δx`

and `Δy`

play a similar role for horizontal grid spacings.

**Example**

```
julia> using Oceananigans
julia> using Oceananigans.AbstractOperations: Δz
julia> c = CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3)));
julia> c_dz = c * Δz # returns BinaryOperation between Field and GridMetricOperation
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
* at (Center, Center, Center)
├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
└── Δzᶜᶜᶜ at (Center, Center, Center)
julia> c .= 1;
julia> c_dz[1, 1, 1]
3.0
```

`Oceananigans.AbstractOperations.Average`

— Method`Average(field::AbstractField; condition = nothing, mask = 0, dims=:)`

Return `Reduction`

representing a spatial average of `field`

over `dims`

.

Over regularly-spaced dimensions this is equivalent to a numerical `mean!`

.

Over dimensions of variable spacing, `field`

is multiplied by the appropriate grid length, area or volume, and divided by the total spatial extent of the interval.

`Oceananigans.AbstractOperations.BinaryOperation`

— Method`BinaryOperation{LX, LY, LZ}(op, a, b, ▶a, ▶b, grid)`

Return an abstract representation of the binary operation `op(▶a(a), ▶b(b))`

on `grid`

, where `▶a`

and `▶b`

interpolate `a`

and `b`

to locations `(LX, LY, LZ)`

.

`Oceananigans.AbstractOperations.ConditionalOperation`

— Method```
ConditionalOperation(operand::AbstractField;
func = identity,
condition = nothing,
mask = 0)
```

Return an abstract representation of a masking procedure applied when `condition`

is satisfied on a field described by `func(operand)`

.

**Positional arguments**

`operand`

: The`AbstractField`

to be masked (it must have a`grid`

property!)

**Keyword arguments**

`func`

: A unary transformation applied element-wise to the field`operand`

at locations where`condition == true`

. Default is`identity`

.`condition`

: either a function of`(i, j, k, grid, operand)`

returning a Boolean, or a 3-dimensional Boolean`AbstractArray`

. At locations where`condition == false`

, operand will be masked by`mask`

`mask`

: the scalar mask

`condition_operand`

is a convenience function used to construct a `ConditionalOperation`

`condition_operand(func::Function, operand::AbstractField, condition, mask) = ConditionalOperation(operand; func, condition, mask)`

**Example**

```
julia> using Oceananigans
julia> using Oceananigans.Fields: condition_operand
julia> c = CenterField(RectilinearGrid(size=(2, 1, 1), extent=(1, 1, 1)));
julia> f(i, j, k, grid, c) = i < 2; d = condition_operand(cos, c, f, 10)
ConditionalOperation at (Center, Center, Center)
├── operand: 2×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── func: cos (generic function with 31 methods)
├── condition: f (generic function with 1 method)
└── mask: 10
julia> d[1, 1, 1]
1.0
julia> d[2, 1, 1]
10
```

`Oceananigans.AbstractOperations.Derivative`

— Method`Derivative{LX, LY, LZ}(∂, arg, ▶, grid)`

Return an abstract representation of the derivative `∂`

on `arg`

, and subsequent interpolation by `▶`

on `grid`

.

`Oceananigans.AbstractOperations.Integral`

— Method`Integral(field::AbstractField; condition = nothing, mask = 0, dims=:)`

Return a `Reduction`

representing a spatial integral of `field`

over `dims`

.

**Example**

Compute the integral of $f(x, y, z) = x y z$ over the domain $(x, y, z) ∈ [0, 1] × [0, 1] × [0, 1]$. The analytical answer is $∭ x y z \, dx \, dy \, dz = 1/8$.

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(8, 8, 8), x=(0, 1), y=(0, 1), z=(0, 1));
julia> f = CenterField(grid);
julia> set!(f, (x, y, z) -> x * y * z)
8×8×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 14×14×14 OffsetArray(::Array{Float64, 3}, -2:11, -2:11, -2:11) with eltype Float64 with indices -2:11×-2:11×-2:11
└── max=0.823975, min=0.000244141, mean=0.125
julia> ∫f = Integral(f)
sum! over dims (1, 2, 3) of BinaryOperation at (Center, Center, Center)
└── operand: BinaryOperation at (Center, Center, Center)
└── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
julia> ∫f = Field(Integral(f));
julia> compute!(∫f);
julia> ∫f[1, 1, 1]
0.125
```

`Oceananigans.AbstractOperations.KernelFunctionOperation`

— Method`KernelFunctionOperation{LX, LY, LZ}(kernel_function, grid, arguments...)`

Construct a `KernelFunctionOperation`

at location `(LX, LY, LZ)`

on `grid`

with `arguments`

.

`kernel_function`

is called with

`kernel_function(i, j, k, grid, arguments...)`

Note that `compute!(kfo::KernelFunctionOperation)`

calls `compute!`

on all `kfo.arguments`

.

**Examples**

Construct a `KernelFunctionOperation`

that returns random numbers:

```
using Oceananigans
grid = RectilinearGrid(size=(1, 8, 8), extent=(1, 1, 1));
random_kernel_function(i, j, k, grid) = rand(); # use CUDA.rand on the GPU
kernel_op = KernelFunctionOperation{Center, Center, Center}(random_kernel_function, grid)
# output
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: random_kernel_function (generic function with 1 method)
└── arguments: ()
```

Construct a `KernelFunctionOperation`

using the vertical vorticity operator used internally to compute vertical vorticity on all grids:

```
using Oceananigans.Operators: ζ₃ᶠᶠᶜ # called with signature ζ₃ᶠᶠᶜ(i, j, k, grid, u, v)
model = HydrostaticFreeSurfaceModel(; grid);
u, v, w = model.velocities;
ζ_op = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, grid, u, v)
# output
KernelFunctionOperation at (Face, Face, Center)
├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: ζ₃ᶠᶠᶜ (generic function with 1 method)
└── arguments: ("1×8×8 Field{Face, Center, Center} on RectilinearGrid on CPU", "1×8×8 Field{Center, Face, Center} on RectilinearGrid on CPU")
```

`Oceananigans.AbstractOperations.UnaryOperation`

— Method`UnaryOperation{LX, LY, LZ}(op, arg, ▶, grid)`

Returns an abstract `UnaryOperation`

representing the action of `op`

on `arg`

, and subsequent interpolation by `▶`

on `grid`

.

`Oceananigans.AbstractOperations.∂x`

— MethodReturn the $x$-derivative function acting at (`X`

, `Y`

, `Any`

).

`Oceananigans.AbstractOperations.∂x`

— Method`∂x(L::Tuple, arg::AbstractField)`

Return an abstract representation of an $x$-derivative acting on field `arg`

followed by interpolation to `L`

, where `L`

is a 3-tuple of `Face`

s and `Center`

s.

`Oceananigans.AbstractOperations.∂x`

— Method`∂x(arg::AbstractField)`

Return an abstract representation of a $x$-derivative acting on field `arg`

.

`Oceananigans.AbstractOperations.∂y`

— MethodReturn the $y$-derivative function acting at (`X`

, `Y`

, `Any`

).

`Oceananigans.AbstractOperations.∂y`

— Method`∂y(L::Tuple, arg::AbstractField)`

Return an abstract representation of a $y$-derivative acting on field `arg`

followed by interpolation to `L`

, where `L`

is a 3-tuple of `Face`

s and `Center`

s.

`Oceananigans.AbstractOperations.∂y`

— Method`∂y(arg::AbstractField)`

Return an abstract representation of a $y$-derivative acting on field `arg`

.

`Oceananigans.AbstractOperations.∂z`

— MethodReturn the $z$-derivative function acting at (`Any`

, `Any`

, `Z`

).

`Oceananigans.AbstractOperations.∂z`

— Method`∂z(L::Tuple, arg::AbstractField)`

Return an abstract representation of a $z$-derivative acting on field `arg`

followed by interpolation to `L`

, where `L`

is a 3-tuple of `Face`

s and `Center`

s.

`Oceananigans.AbstractOperations.∂z`

— Method`∂z(arg::AbstractField)`

Return an abstract representation of a $z$-derivative acting on field `arg`

.

`Oceananigans.AbstractOperations.@at`

— Macro`@at location abstract_operation`

Modify the `abstract_operation`

so that it returns values at `location`

, where `location`

is a 3-tuple of `Face`

s and `Center`

s.

`Oceananigans.AbstractOperations.@binary`

— Macro`@binary op1 op2 op3...`

Turn each binary function in the list `(op1, op2, op3...)`

into a binary operator on `Oceananigans.Fields`

for use in `AbstractOperations`

.

Note: a binary function is a function with two arguments: for example, `+(x, y)`

is a binary function.

Also note: a binary function in `Base`

must be imported to be extended: use `import Base: op; @binary op`

.

**Example**

```
julia> using Oceananigans, Oceananigans.AbstractOperations
julia> using Oceananigans.AbstractOperations: BinaryOperation, AbstractGridMetric, choose_location
julia> plus_or_times(x, y) = x < 0 ? x + y : x * y
plus_or_times (generic function with 1 method)
julia> @binary plus_or_times
Set{Any} with 6 elements:
:+
:/
:^
:-
:*
:plus_or_times
julia> c, d = (CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) for i = 1:2);
julia> plus_or_times(c, d)
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
plus_or_times at (Center, Center, Center)
├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
└── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
```

`Oceananigans.AbstractOperations.@multiary`

— Macro`@multiary op1 op2 op3...`

Turn each multiary operator in the list `(op1, op2, op3...)`

into a multiary operator on `Oceananigans.Fields`

for use in `AbstractOperations`

.

Note that a multiary operator:

- is a function with two or more arguments: for example,
`+(x, y, z)`

is a multiary function; - must be imported to be extended if part of
`Base`

: use`import Base: op; @multiary op`

; - can only be called on
`Oceananigans.Field`

s if the "location" is noted explicitly; see example.

**Example**

```
julia> using Oceananigans, Oceananigans.AbstractOperations
julia> harmonic_plus(a, b, c) = 1/3 * (1/a + 1/b + 1/c)
harmonic_plus (generic function with 1 method)
julia> c, d, e = Tuple(CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) for i = 1:3);
julia> harmonic_plus(c, d, e) # before magic @multiary transformation
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
* at (Center, Center, Center)
├── 0.3333333333333333
└── + at (Center, Center, Center)
├── / at (Center, Center, Center)
│ ├── 1
│ └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── / at (Center, Center, Center)
│ ├── 1
│ └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
└── / at (Center, Center, Center)
├── 1
└── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
julia> @multiary harmonic_plus
Set{Any} with 3 elements:
:+
:harmonic_plus
:*
julia> harmonic_plus(c, d, e)
MultiaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
harmonic_plus at (Center, Center, Center)
├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
└── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
```

`Oceananigans.AbstractOperations.@unary`

— Macro`@unary op1 op2 op3...`

Turn each unary function in the list `(op1, op2, op3...)`

into a unary operator on `Oceananigans.Fields`

for use in `AbstractOperations`

.

Note: a unary function is a function with one argument: for example, `sin(x)`

is a unary function.

Also note: a unary function in `Base`

must be imported to be extended: use `import Base: op; @unary op`

.

**Example**

```
julia> using Oceananigans, Oceananigans.Grids, Oceananigans.AbstractOperations
julia> square_it(x) = x^2
square_it (generic function with 1 method)
julia> @unary square_it
Set{Any} with 10 elements:
:+
:sqrt
:square_it
:cos
:exp
:interpolate_identity
:-
:tanh
:sin
:abs
julia> c = CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)));
julia> square_it(c)
UnaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
square_it at (Center, Center, Center) via identity
└── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
```

## Advection

`Oceananigans.Advection.Centered`

— Type`struct Centered{N, FT, XT, YT, ZT, CA} <: AbstractCenteredAdvectionScheme{N, FT}`

Centered reconstruction scheme.

`Oceananigans.Advection.UpwindBiased`

— Type`struct UpwindBiasedFifthOrder <: AbstractUpwindBiasedAdvectionScheme{3}`

Upwind-biased fifth-order advection scheme.

`Oceananigans.Advection.VectorInvariant`

— Method```
VectorInvariant(; vorticity_scheme = EnstrophyConserving(),
vorticity_stencil = VelocityStencil(),
vertical_scheme = EnergyConserving(),
kinetic_energy_gradient_scheme = vertical_scheme)
```

Return a vector invariant momentum advection scheme.

**Keyword arguments**

`vorticity_scheme`

: Scheme used for`Center`

reconstruction of vorticity. Default:`EnstrophyConserving()`

. Options: -`UpwindBiased()`

-`WENO()`

-`EnergyConserving()`

-`EnstrophyConserving()`

`vorticity_stencil`

: Stencil used for smoothness indicators for`WENO`

schemes. Default:`VelocityStencil()`

. Options: -`VelocityStencil()`

(smoothness based on horizontal velocities) -`DefaultStencil()`

(smoothness based on variable being reconstructed)`vertical_scheme`

: Scheme used for vertical advection of horizontal momentum. Default:`EnergyConserving()`

.`kinetic_energy_gradient_scheme`

: Scheme used for kinetic energy gradient reconstruction. Default:`vertical_scheme`

.`divergence_scheme`

: Scheme used for divergence flux. Only upwinding schemes are supported. Default:`vorticity_scheme`

.`upwinding`

: Treatment of upwinded reconstruction of divergence and kinetic energy gradient. Default:`OnlySelfUpwinding()`

. Options: -`CrossAndSelfUpwinding()`

-`OnlySelfUpwinding()`

-`VelocityUpwinding()`

`multi_dimensional_stencil`

: whether or not to use a horizontal two-dimensional stencil for the reconstruction of vorticity, divergence and kinetic energy gradient. Currently the "tangential" direction uses 5th-order centered WENO reconstruction.

**Examples**

```
julia> using Oceananigans
julia> VectorInvariant()
Vector Invariant, Dimension-by-dimension reconstruction
Vorticity flux scheme:
└── EnstrophyConserving{Float64}
Vertical advection / Divergence flux scheme:
└── EnergyConserving{Float64}
```

```
julia> using Oceananigans
julia> VectorInvariant(vorticity_scheme = WENO(), vertical_scheme = WENO(order = 3))
Vector Invariant, Dimension-by-dimension reconstruction
Vorticity flux scheme:
├── WENO reconstruction order 5
└── smoothness ζ: Oceananigans.Advection.VelocityStencil()
Vertical advection / Divergence flux scheme:
├── WENO reconstruction order 3
└── upwinding treatment: OnlySelfUpwinding
KE gradient and Divergence flux cross terms reconstruction:
└── Centered reconstruction order 2
Smoothness measures:
├── smoothness δU: FunctionStencil f = divergence_smoothness
├── smoothness δV: FunctionStencil f = divergence_smoothness
├── smoothness δu²: FunctionStencil f = u_smoothness
└── smoothness δv²: FunctionStencil f = v_smoothness
```

`Oceananigans.Advection.WENO`

— Type```
WENO([FT=Float64;]
order = 5,
grid = nothing,
zweno = true,
bounds = nothing)
```

Construct a weigthed essentially non-oscillatory advection scheme of order `order`

.

**Keyword arguments**

`order`

: The order of the WENO advection scheme. Default: 5`grid`

: (defaults to`nothing`

)`zweno`

: When`true`

implement a Z-WENO formulation for the WENO weights calculation. (defaults to`true`

)

**Examples**

```
julia> using Oceananigans
julia> WENO()
WENO reconstruction order 5
Smoothness formulation:
└── Z-weno
Boundary scheme:
└── WENO reconstruction order 3
Symmetric scheme:
└── Centered reconstruction order 4
Directions:
├── X regular
├── Y regular
└── Z regular
```

```
julia> using Oceananigans
julia> Nx, Nz = 16, 10;
julia> Lx, Lz = 1e4, 1e3;
julia> chebychev_spaced_z_faces(k) = - Lz/2 - Lz/2 * cos(π * (k - 1) / Nz);
julia> grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4), topology=(Periodic, Flat, Bounded),
x = (0, Lx), z = chebychev_spaced_z_faces);
julia> WENO(grid; order=7)
WENO reconstruction order 7
Smoothness formulation:
└── Z-weno
Boundary scheme:
└── WENO reconstruction order 5
Symmetric scheme:
└── Centered reconstruction order 6
Directions:
├── X regular
├── Y regular
└── Z stretched
```

`Oceananigans.Advection.WENOVectorInvariant`

— Method```
function WENOVectorInvariant(; upwinding = nothing,
multi_dimensional_stencil = false,
weno_kw...)
```

`Oceananigans.Advection.div_Uc`

— Method`div_uc(i, j, k, grid, advection, U, c)`

Calculate the divergence of the flux of a tracer quantity $c$ being advected by a velocity field, $𝛁⋅(𝐯 c)$,

`1/V * [δxᶜᵃᵃ(Ax * u * ℑxᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * v * ℑyᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * w * ℑzᵃᵃᶠ(c))]`

which ends up at the location `ccc`

.

`Oceananigans.Advection.div_𝐯u`

— Method`div_𝐯u(i, j, k, grid, advection, U, u)`

Calculate the advection of momentum in the $x$-direction using the conservative form, $𝛁⋅(𝐯 u)$,

`1/Vᵘ * [δxᶠᵃᵃ(ℑxᶜᵃᵃ(Ax * u) * ℑxᶜᵃᵃ(u)) + δy_fca(ℑxᶠᵃᵃ(Ay * v) * ℑyᵃᶠᵃ(u)) + δz_fac(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(u))]`

which ends up at the location `fcc`

.

`Oceananigans.Advection.div_𝐯v`

— Method`div_𝐯v(i, j, k, grid, advection, U, v)`

Calculate the advection of momentum in the $y$-direction using the conservative form, $𝛁⋅(𝐯 v)$,

`1/Vʸ * [δx_cfa(ℑyᵃᶠᵃ(Ax * u) * ℑxᶠᵃᵃ(v)) + δyᵃᶠᵃ(ℑyᵃᶜᵃ(Ay * v) * ℑyᵃᶜᵃ(v)) + δz_afc(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(w))]`

which ends up at the location `cfc`

.

`Oceananigans.Advection.div_𝐯w`

— Method`div_𝐯w(i, j, k, grid, advection, U, w)`

Calculate the advection of momentum in the $z$-direction using the conservative form, $𝛁⋅(𝐯 w)$,

`1/Vʷ * [δx_caf(ℑzᵃᵃᶠ(Ax * u) * ℑxᶠᵃᵃ(w)) + δy_acf(ℑzᵃᵃᶠ(Ay * v) * ℑyᵃᶠᵃ(w)) + δzᵃᵃᶠ(ℑzᵃᵃᶜ(Az * w) * ℑzᵃᵃᶜ(w))]`

which ends up at the location `ccf`

.

## Architectures

`Oceananigans.Architectures.AbstractArchitecture`

— Type`AbstractArchitecture`

Abstract supertype for architectures supported by Oceananigans.

`Oceananigans.Architectures.CPU`

— Type`CPU <: AbstractArchitecture`

Run Oceananigans on one CPU node. Uses multiple threads if the environment variable `JULIA_NUM_THREADS`

is set.

`Oceananigans.Architectures.GPU`

— Type`GPU <: AbstractArchitecture`

Run Oceananigans on a single NVIDIA CUDA GPU.

## Boundary conditions

`Oceananigans.BoundaryConditions.BoundaryCondition`

— Type`struct BoundaryCondition{C<:AbstractBoundaryConditionClassification, T}`

Container for boundary conditions.

`Oceananigans.BoundaryConditions.BoundaryCondition`

— Method`BoundaryCondition(Classification::DataType, condition)`

Construct a boundary condition of type `BC`

with a number or array as a `condition`

.

Boundary condition types include `Periodic`

, `Flux`

, `Value`

, `Gradient`

, and `Open`

.

`Oceananigans.BoundaryConditions.BoundaryCondition`

— Method```
BoundaryCondition(Classification::DataType, condition::Function;
parameters = nothing,
discrete_form = false,
field_dependencies=())
```

Construct a boundary condition of type `Classification`

with a function boundary `condition`

.

By default, the function boudnary `condition`

is assumed to have the 'continuous form' `condition(ξ, η, t)`

, where `t`

is time and `ξ`

and `η`

vary along the boundary. In particular:

- On
`x`

-boundaries,`condition(y, z, t)`

. - On
`y`

-boundaries,`condition(x, z, t)`

. - On
`z`

-boundaries,`condition(x, y, t)`

.

If `parameters`

is not `nothing`

, then function boundary conditions have the form `func(ξ, η, t, parameters)`

, where `ξ`

and `η`

are spatial coordinates varying along the boundary as explained above.

If `discrete_form = true`

, the function `condition`

is assumed to have the "discrete form",

`condition(i, j, grid, clock, model_fields)`

where `i`

, and `j`

are indices that vary along the boundary. If `discrete_form = true`

and `parameters`

is not `nothing`

, the function `condition`

is called with

`condition(i, j, grid, clock, model_fields, parameters)`

`Oceananigans.BoundaryConditions.FieldBoundaryConditions`

— Type```
FieldBoundaryConditions(grid, location, indices=(:, :, :);
west = default_auxiliary_bc(topology(grid, 1)(), location[1]()),
east = default_auxiliary_bc(topology(grid, 1)(), location[1]()),
south = default_auxiliary_bc(topology(grid, 2)(), location[2]()),
north = default_auxiliary_bc(topology(grid, 2)(), location[2]()),
bottom = default_auxiliary_bc(topology(grid, 3)(), location[3]()),
top = default_auxiliary_bc(topology(grid, 3)(), location[3]()),
immersed = NoFluxBoundaryCondition())
```

Return boundary conditions for auxiliary fields (fields whose values are derived from a model's prognostic fields) on `grid`

and at `location`

.

**Keyword arguments**

Keyword arguments specify boundary conditions on the 6 possible boundaries:

`west`

, left end point in the`x`

-direction where`i = 1`

`east`

, right end point in the`x`

-direction where`i = grid.Nx`

`south`

, left end point in the`y`

-direction where`j = 1`

`north`

, right end point in the`y`

-direction where`j = grid.Ny`

`bottom`

, right end point in the`z`

-direction where`k = 1`

`top`

, right end point in the`z`

-direction where`k = grid.Nz`

`immersed`

: boundary between solid and fluid for immersed boundaries

If a boundary condition is unspecified, the default for auxiliary fields and the topology in the boundary-normal direction is used:

`PeriodicBoundaryCondition`

for`Periodic`

directions`GradientBoundaryCondition(0)`

for`Bounded`

directions and`Centered`

-located fields`nothing`

for`Bounded`

directions and`Face`

-located fields`nothing`

for`Flat`

directions and/or`Nothing`

-located fields

`Oceananigans.BoundaryConditions.FieldBoundaryConditions`

— Type`FieldBoundaryConditions(; kwargs...)`

Return a template for boundary conditions on prognostic fields.

**Keyword arguments**

Keyword arguments specify boundary conditions on the 7 possible boundaries:

`west`

: left end point in the`x`

-direction where`i = 1`

`east`

: right end point in the`x`

-direction where`i = grid.Nx`

`south`

: left end point in the`y`

-direction where`j = 1`

`north`

: right end point in the`y`

-direction where`j = grid.Ny`

`bottom`

: right end point in the`z`

-direction where`k = 1`

`top`

: right end point in the`z`

-direction where`k = grid.Nz`

`immersed`

: boundary between solid and fluid for immersed boundaries

If a boundary condition is unspecified, the default for prognostic fields and the topology in the boundary-normal direction is used:

`PeriodicBoundaryCondition`

for`Periodic`

directions`NoFluxBoundaryCondition`

for`Bounded`

directions and`Centered`

-located fields`ImpenetrableBoundaryCondition`

for`Bounded`

directions and`Face`

-located fields`nothing`

for`Flat`

directions and/or`Nothing`

-located fields

`Oceananigans.BoundaryConditions.Flux`

— Type`struct Flux <: AbstractBoundaryConditionClassification`

A classification specifying a boundary condition on the flux of a field.

The sign convention is such that a positive flux represents the flux of a quantity in the positive direction. For example, a positive vertical flux implies a quantity is fluxed upwards, in the $+z$ direction.

Due to this convention, a positive flux applied to the top boundary specifies that a quantity is fluxed upwards across the top boundary and thus out of the domain. As a result, a positive flux applied to a top boundary leads to a reduction of that quantity in the interior of the domain; for example, a positive, upwards flux of heat at the top of the domain acts to cool the interior of the domain. Conversely, a positive flux applied to the bottom boundary leads to an increase of the quantity in the interior of the domain. The same logic holds for east, west, north, and south boundaries.

`Oceananigans.BoundaryConditions.Gradient`

— Type`struct Gradient <: AbstractBoundaryConditionClassification`

A classification specifying a boundary condition on the derivative or gradient of a field. Also called a Neumann boundary condition.

`Oceananigans.BoundaryConditions.Open`

— Type`struct Open <: AbstractBoundaryConditionClassification`

A classification that specifies the halo regions of a field directly.

For fields located at Faces, Open also specifies field value *on* the boundary.

Open boundary conditions are used to specify the component of a velocity field normal to a boundary and can also be used to describe nested or linked simulation domains.

`Oceananigans.BoundaryConditions.Value`

— Type`struct Value <: AbstractBoundaryConditionClassification`

A classification specifying a boundary condition on the value of a field. Also called a Dirchlet boundary condition.

`Oceananigans.BoundaryConditions.apply_x_bcs!`

— MethodApply flux boundary conditions to a field `c`

by adding the associated flux divergence to the source term `Gc`

at the left and right.

`Oceananigans.BoundaryConditions.apply_y_bcs!`

— MethodApply flux boundary conditions to a field `c`

by adding the associated flux divergence to the source term `Gc`

at the left and right.

`Oceananigans.BoundaryConditions.apply_z_bcs!`

— MethodApply flux boundary conditions to a field `c`

by adding the associated flux divergence to the source term `Gc`

at the top and bottom.

`Oceananigans.BoundaryConditions.fill_halo_regions!`

— MethodFill halo regions in $x$, $y$, and $z$ for a given field's data.

## Buoyancy models

`Oceananigans.BuoyancyModels.Buoyancy`

— Method`Buoyancy(; model, gravity_unit_vector=NegativeZDirection())`

Construct a `buoyancy`

given a buoyancy `model`

. Optional keyword argument `gravity_unit_vector`

can be used to specify the direction of gravity (default `NegativeZDirection()`

). The buoyancy acceleration acts in the direction opposite to gravity.

**Example**

```
using Oceananigans
grid = RectilinearGrid(size=(1, 8, 8), extent=(1, 1, 1))
θ = 45 # degrees
g̃ = (0, -sind(θ), -cosd(θ))
buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃)
model = NonhydrostaticModel(; grid, buoyancy, tracers=:b)
# output
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = Tuple{Float64, Float64, Float64}
└── coriolis: Nothing
```

`Oceananigans.BuoyancyModels.BuoyancyTracer`

— Type`BuoyancyTracer <: AbstractBuoyancyModel{Nothing}`

Type indicating that the tracer `b`

represents buoyancy.

`Oceananigans.BuoyancyModels.LinearEquationOfState`

— Type`LinearEquationOfState([FT=Float64;] thermal_expansion=1.67e-4, haline_contraction=7.80e-4)`

Return `LinearEquationOfState`

for `SeawaterBuoyancy`

with `thermal_expansion`

coefficient and `haline_contraction`

coefficient. The buoyancy perturbation $b$ for `LinearEquationOfState`

is

\[ b = g (α T - β S),\]

where $g$ is gravitational acceleration, $α$ is `thermal_expansion`

, $β$ is `haline_contraction`

, $T$ is temperature, and $S$ is practical salinity units.

Default constants in units inverse Kelvin and practical salinity units for `thermal_expansion`

and `haline_contraction`

, respectively, are taken from Table 1.2 (page 33) of Vallis, "Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation" (2nd ed, 2017).

`Oceananigans.BuoyancyModels.LinearEquationOfState`

— Type`LinearEquationOfState{FT} <: AbstractEquationOfState`

Linear equation of state for seawater.

`Oceananigans.BuoyancyModels.SeawaterBuoyancy`

— Type```
SeawaterBuoyancy([FT = Float64;]
gravitational_acceleration = g_Earth,
equation_of_state = LinearEquationOfState(FT),
constant_temperature = nothing,
constant_salinity = nothing)
```

Return parameters for a temperature- and salt-stratified seawater buoyancy model with a `gravitational_acceleration`

constant (typically called $g$), and an `equation_of_state`

that related temperature and salinity (or conservative temperature and absolute salinity) to density anomalies and buoyancy.

Setting `constant_temperature`

to something that is not `nothing`

indicates that buoyancy depends only on salinity. For a nonlinear equation of state, the value provided `constant_temperature`

is used as the temperature of the system. Vice versa, setting `constant_salinity`

indicates that buoyancy depends only on temperature.

For a linear equation of state, the values of `constant_temperature`

or `constant_salinity`

are irrelevant.

**Examples**

The "TEOS10" equation of state, see https://www.teos-10.org

```
julia> using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
julia> teos10 = TEOS10EquationOfState()
BoussinesqEquationOfState{Float64}:
├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64}
└── reference_density: 1020.0
```

Buoyancy that depends on both temperature and salinity

```
julia> using Oceananigans
julia> buoyancy = SeawaterBuoyancy(equation_of_state=teos10)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}
```

Buoyancy that depends only on salinity with temperature held at 20 degrees Celsius

```
julia> salinity_dependent_buoyancy = SeawaterBuoyancy(equation_of_state=teos10, constant_temperature=20)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
├── constant_temperature: 20
└── equation_of_state: BoussinesqEquationOfState{Float64}
```

Buoyancy that depends only on temperature with salinity held at 35 psu

```
julia> temperature_dependent_buoyancy = SeawaterBuoyancy(equation_of_state=teos10, constant_salinity=35)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
├── constant_salinity: 35
└── equation_of_state: BoussinesqEquationOfState{Float64}
```

`Oceananigans.BuoyancyModels.SeawaterBuoyancy`

— Type`SeawaterBuoyancy{FT, EOS, T, S} <: AbstractBuoyancyModel{EOS}`

BuoyancyModels model for seawater. `T`

and `S`

are either `nothing`

if both temperature and salinity are active, or of type `FT`

if temperature or salinity are constant, respectively.

`Oceananigans.BuoyancyModels.∂x_b`

— Method`∂x_b(i, j, k, grid, b::SeawaterBuoyancy, C)`

Returns the $x$-derivative of buoyancy for temperature and salt-stratified water,

\[∂_x b = g ( α ∂_x T - β ∂_x S ) ,\]

where $g$ is gravitational acceleration, $α$ is the thermal expansion coefficient, $β$ is the haline contraction coefficient, $T$ is conservative temperature, and $S$ is absolute salinity.

Note: In Oceananigans, `model.tracers.T`

is conservative temperature and `model.tracers.S`

is absolute salinity.

Note that $∂_x T$ (`∂x_T`

), $∂_x S$ (`∂x_S`

), $α$, and $β$ are all evaluated at cell interfaces in `x`

and cell centers in `y`

and `z`

.

`Oceananigans.BuoyancyModels.∂y_b`

— Method`∂y_b(i, j, k, grid, b::SeawaterBuoyancy, C)`

Returns the $y$-derivative of buoyancy for temperature and salt-stratified water,

\[∂_y b = g ( α ∂_y T - β ∂_y S ) ,\]

where $g$ is gravitational acceleration, $α$ is the thermal expansion coefficient, $β$ is the haline contraction coefficient, $T$ is conservative temperature, and $S$ is absolute salinity.

Note: In Oceananigans, `model.tracers.T`

is conservative temperature and `model.tracers.S`

is absolute salinity.

Note that $∂_y T$ (`∂y_T`

), $∂_y S$ (`∂y_S`

), $α$, and $β$ are all evaluated at cell interfaces in `y`

and cell centers in `x`

and `z`

.

`Oceananigans.BuoyancyModels.∂z_b`

— Method`∂z_b(i, j, k, grid, b::SeawaterBuoyancy, C)`

Returns the vertical derivative of buoyancy for temperature and salt-stratified water,

\[∂_z b = N^2 = g ( α ∂_z T - β ∂_z S ) ,\]

where $g$ is gravitational acceleration, $α$ is the thermal expansion coefficient, $β$ is the haline contraction coefficient, $T$ is conservative temperature, and $S$ is absolute salinity.

Note: In Oceananigans, `model.tracers.T`

is conservative temperature and `model.tracers.S`

is absolute salinity.

Note that $∂_z T$ (`∂z_T`

), $∂_z S$ (`∂z_S`

), $α$, and $β$ are all evaluated at cell interfaces in `z`

and cell centers in `x`

and `y`

.

## Coriolis

`Oceananigans.Coriolis.ActiveCellEnstrophyConserving`

— Type`struct ActiveCellEnstrophyConserving`

A parameter object for an enstrophy-conserving Coriolis scheme that excludes inactive (dry/land) edges (indices for which `peripheral_node == true`

) from the velocity interpolation.

`Oceananigans.Coriolis.BetaPlane`

— Type```
BetaPlane([T=Float64;] f₀=nothing, β=nothing,
rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)
```

Return a $β$-plane Coriolis parameter, $f = f₀ + β y$.

The user may specify both `f₀`

and `β`

, or the three parameters `rotation_rate`

, `latitude`

(in degrees), and `radius`

that specify the rotation rate and radius of a planet, and the central latitude (where $y = 0$) at which the `β`

-plane approximation is to be made.

If `f₀`

and `β`

are not specified, they are calculated from `rotation_rate`

, `latitude`

, and `radius`

according to the relations `f₀ = 2 * rotation_rate * sind(latitude)`

and `β = 2 * rotation_rate * cosd(latitude) / radius`

.

By default, the `rotation_rate`

and planet `radius`

are assumed to be Earth's.

`Oceananigans.Coriolis.BetaPlane`

— Type`struct BetaPlane{T} <: AbstractRotation`

A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + β y`

) that accounts for the variation of the locally vertical component of the rotation vector with latitude.

`Oceananigans.Coriolis.ConstantCartesianCoriolis`

— Type```
ConstantCartesianCoriolis([FT=Float64;] fx=nothing, fy=nothing, fz=nothing,
f=nothing, rotation_axis=ZDirection(),
rotation_rate=Ω_Earth, latitude=nothing)
```

Return a parameter object for a constant rotation decomposed into the `x`

, `y`

, and `z`

directions. In oceanography the components `x`

, `y`

, `z`

correspond to the directions east, north, and up. This constant rotation can be specified in three different ways:

- Specifying all components
`fx`

,`fy`

and`fz`

directly. - Specifying the Coriolis parameter
`f`

and (optionally) a`rotation_axis`

(which defaults to the`z`

direction if not specified). - Specifying
`latitude`

(in degrees) and (optionally) a`rotation_rate`

in radians per second (which defaults to Earth's rotation rate).

`Oceananigans.Coriolis.ConstantCartesianCoriolis`

— Type`struct ConstantCartesianCoriolis{FT} <: AbstractRotation`

A Coriolis implementation that accounts for the locally vertical and possibly both local horizontal components of a constant rotation vector. This is a more general implementation of `FPlane`

, which only accounts for the locally vertical component.

`Oceananigans.Coriolis.FPlane`

— Type`FPlane([FT=Float64;] f=nothing, rotation_rate=Ω_Earth, latitude=nothing)`

Return a parameter object for constant rotation at the angular frequency `f/2`

, and therefore with background vorticity `f`

, around a vertical axis. If `f`

is not specified, it is calculated from `rotation_rate`

and `latitude`

(in degrees) according to the relation `f = 2 * rotation_rate * sind(latitude)`

.

By default, `rotation_rate`

is assumed to be Earth's.

Also called `FPlane`

, after the "f-plane" approximation for the local effect of a planet's rotation in a planar coordinate system tangent to the planet's surface.

`Oceananigans.Coriolis.FPlane`

— Type`struct FPlane{FT} <: AbstractRotation`

A parameter object for constant rotation around a vertical axis.

`Oceananigans.Coriolis.HydrostaticSphericalCoriolis`

— Type`struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation`

A parameter object for constant rotation around a vertical axis on the sphere.

`Oceananigans.Coriolis.HydrostaticSphericalCoriolis`

— Method```
HydrostaticSphericalCoriolis([FT=Float64;]
rotation_rate = Ω_Earth,
scheme = EnergyConserving())
```

Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`

. By default, `rotation_rate`

is assumed to be Earth's.

**Keyword arguments**

`scheme`

: Either`EnergyConserving()`

(default),`EnstrophyConserving()`

, or`ActiveCellEnstrophyConserving()`

.

`Oceananigans.Coriolis.NonTraditionalBetaPlane`

— Type```
NonTraditionalBetaPlane(FT=Float64;
fz=nothing, fy=nothing, β=nothing, γ=nothing,
rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)
```

The user may directly specify `fz`

, `fy`

, `β`

, `γ`

, and `radius`

or the three parameters `rotation_rate`

, `latitude`

(in degrees), and `radius`

that specify the rotation rate and radius of a planet, and the central latitude (where $y = 0$) at which the non-traditional `β`

-plane approximation is to be made.

If `fz`

, `fy`

, `β`

, and `γ`

are not specified, they are calculated from `rotation_rate`

, `latitude`

, and `radius`

according to the relations `fz = 2 * rotation_rate * sind(latitude)`

, `fy = 2 * rotation_rate * cosd(latitude)`

, `β = 2 * rotation_rate * cosd(latitude) / radius`

, and `γ = - 4 * rotation_rate * sind(latitude) / radius`

.

By default, the `rotation_rate`

and planet `radius`

is assumed to be Earth's.

`Oceananigans.Coriolis.NonTraditionalBetaPlane`

— Type`struct NonTraditionalBetaPlane{FT} <: AbstractRotation`

A Coriolis implementation that accounts for the latitudinal variation of both the locally vertical and the locally horizontal components of the rotation vector. The "traditional" approximation in ocean models accounts for only the locally vertical component of the rotation vector (see `BetaPlane`

).

This implementation is based off of section 5 of Dellar (2011) and it conserves energy, angular momentum, and potential vorticity.

**References**

Dellar, P. (2011). Variations on a beta-plane: Derivation of non-traditional beta-plane equations from Hamilton's principle on a sphere. Journal of Fluid Mechanics, 674, 174-195. doi:10.1017/S0022112010006464

## Diagnostics

`Oceananigans.Diagnostics.CFL`

— Type`struct CFL{D, S}`

An object for computing the Courant-Freidrichs-Lewy (CFL) number.

`Oceananigans.Diagnostics.CFL`

— Method`CFL(Δt [, timescale = Oceananigans.Advection.cell_advection_timescale])`

Return an object for computing the Courant-Freidrichs-Lewy (CFL) number associated with time step `Δt`

or `TimeStepWizard`

and `timescale`

.

See also `AdvectiveCFL`

and `DiffusiveCFL`

.

`Oceananigans.Diagnostics.StateChecker`

— Method`StateChecker(; schedule, fields)`

Returns a `StateChecker`

that logs field information (minimum, maximum, mean) for each field in a named tuple of `fields`

when `schedule`

actuates.

`Oceananigans.Diagnostics.AdvectiveCFL`

— Method`AdvectiveCFL(Δt)`

Return an object for computing the Courant-Freidrichs-Lewy (CFL) number associated with time step `Δt`

or `TimeStepWizard`

and the time scale for advection across a cell. The advective CFL is, e.g., $U Δt / Δx$.

**Example**

```
julia> using Oceananigans
julia> model = NonhydrostaticModel(grid = RectilinearGrid(size=(16, 16, 16), extent=(8, 8, 8)));
julia> Δt = 1.0;
julia> cfl = AdvectiveCFL(Δt);
julia> model.velocities.u .= π;
julia> cfl(model)
6.283185307179586
```

`Oceananigans.Diagnostics.DiffusiveCFL`

— Method`DiffusiveCFL(Δt)`

Returns an object for computing the diffusive Courant-Freidrichs-Lewy (CFL) number associated with time step `Δt`

or `TimeStepWizard`

and the time scale for diffusion across a cell associated with `model.closure`

. The diffusive CFL, e.g., for viscosity is $ν Δt / Δx²$.

The maximum diffusive CFL number among viscosity and all tracer diffusivities is returned.

**Example**

```
julia> using Oceananigans
julia> model = NonhydrostaticModel(grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1)),
closure = ScalarDiffusivity(; ν = 1e-2));
julia> Δt = 0.1;
julia> dcfl = DiffusiveCFL(Δt);
julia> dcfl(model)
0.256
```

## Distributed

`Oceananigans.DistributedComputations.Distributed`

— Type```
Distributed(child_architecture = CPU();
topology,
partition,
devices = nothing,
communicator = MPI.COMM_WORLD)
```

Constructor for a distributed architecture that uses MPI for communications

**Positional arguments**

`child_architecture`

: Specifies whether the computation is performed on CPUs or GPUs. Default:`child_architecture = CPU()`

.

**Keyword arguments**

`synchronized_communication`

: if true, always use synchronized communication through ranks`ranks`

(required): A 3-tuple`(Rx, Ry, Rz)`

specifying the total processors in the`x`

,`y`

and`z`

direction. NOTE: support for distributed z direction is limited, so`Rz = 1`

is strongly suggested.`devices`

:`GPU`

device linked to local rank. The GPU will be assigned based on the local node rank as such`devices[node_rank]`

. Make sure to run`--ntasks-per-node`

<=`--gres=gpu`

. If`nothing`

, the devices will be assigned automatically based on the available resources`communicator`

: the MPI communicator,`MPI.COMM_WORLD`

. This keyword argument should not be tampered with if not for testing or developing. Change at your own risk!

`Oceananigans.DistributedComputations.DistributedFFTBasedPoissonSolver`

— Method`DistributedFFTBasedPoissonSolver(global_grid, local_grid)`

Return a FFT-based solver for the Poisson equation,

\[∇²φ = b\]

for `Distributed`

itectures.

**Supported configurations**

We support two "modes":

Vertical pencil decompositions: two-dimensional decompositions in $(x, y)$ for three dimensional problems that satisfy either

`Nz > Rx`

or`Nz > Ry`

.One-dimensional decompositions in either $x$ or $y$.

Above, `Nz = size(global_grid, 3)`

and `Rx, Ry, Rz = architecture(local_grid).ranks`

.

Other configurations that are decomposed in $(x, y)$ but have too few `Nz`

, or any configuration decomposed in $z$, are *not* supported.

**Algorithm for two-dimensional decompositions**

For two-dimensional decompositions (useful for three-dimensional problems), there are three forward transforms, three backward transforms, and four transpositions requiring MPI communication. In the schematic below, the first dimension is always the local dimension. In our implementation of the PencilFFTs algorithm, we require *either* `Nz >= Rx`

, *or* `Nz >= Ry`

, where `Nz`

is the number of vertical cells, `Rx`

is the number of ranks in $x$, and `Ry`

is the number of ranks in $y$. Below, we outline the algorithm for the case `Nz >= Rx`

. If `Nz < Rx`

, but `Nz > Ry`

, a similar algorithm applies with $x$ and $y$ swapped:

`first(storage)`

is initialized with layout $(z, x, y)$.- Transform along $z$.

3 Transpose + communicate to `storage[2]`

in layout $(x, z, y)$, which is distributed into `(Rx, Ry)`

processes in $(z, y)$.

- Transform along $x$.
- Transpose + communicate to
`last(storage)`

in layout $(y, x, z)$, which is distributed into`(Rx, Ry)`

processes in $(x, z)$. - Transform in $y$.

At this point the three in-place forward transforms are complete, and we solve the Poisson equation by updating `last(storage)`

. Then the process is reversed to obtain `first(storage)`

in physical space and with the layout $(z, x, y)$.

**Restrictions**

The algorithm for two-dimensional decompositions requires that `Nz = size(global_grid, 3)`

is larger than either `Rx = ranks[1]`

or `Ry = ranks[2]`

, where `ranks`

are configured when building `Distributed`

. If `Nz`

does not satisfy this condition, we can only support a one-dimensional decomposition.

**Algorithm for one-dimensional decompositions**

This algorithm requires a one-dimensional decomposition with *either* `Rx = 1`

*or* `Ry = 1`

, and is important to support two-dimensional transforms.

For one-dimensional decompositions, we place the decomposed direction *last*. If the number of ranks is `Rh = max(Rx, Ry)`

, this algorithm requires that *both* `Nx > Rh`

*and* `Ny > Rh`

. The resulting flow of transposes and transforms is similar to the two-dimensional case. It remains somewhat of a mystery why this succeeds (i.e., why the last transform is correctly decomposed).

`Oceananigans.DistributedComputations.Equal`

— Type`Equal()`

Return a type that partitions a direction equally among remaining processes.

`Equal()`

can be used for only one direction. Other directions must either be unspecified, or specifically defined by `Int`

, `Fractional`

, or `Sizes`

.

`Oceananigans.DistributedComputations.Fractional`

— Method`Fractional(ϵ₁, ϵ₂, ..., ϵₙ)`

Return a type that partitions a direction unequally. The total work is `W = sum(ϵᵢ)`

, and each process is then allocated `ϵᵢ / W`

of the domain.

`Oceananigans.DistributedComputations.Partition`

— Method`Partition(; x = 1, y = 1, z = 1)`

Return `Partition`

representing the division of a domain in the `x`

(first), `y`

(second) and `z`

(third) dimension

**Keyword arguments:**

`x`

: partitioning of the first dimension`y`

: partitioning of the second dimension`z`

: partitioning of the third dimension

if supplied as positional arguments `x`

will be the first argument, `y`

the second and `z`

the third

`x`

, `y`

and `z`

can be:

`x::Int`

: allocate`x`

processors to the first dimension`Equal()`

: divide the domain in`x`

equally among the remaining processes (not supported for multiple directions)`Fractional(ϵ₁, ϵ₂, ..., ϵₙ):`

divide the domain unequally among`N`

processes. The total work is`W = sum(ϵᵢ)`

, and each process is then allocated`ϵᵢ / W`

of the domain.`Sizes(n₁, n₂, ..., nₙ)`

: divide the domain unequally. The total work is`W = sum(nᵢ)`

, and each process is then allocated`nᵢ`

.

**Examples:**

```
julia> using Oceananigans; using Oceananigans.DistributedComputations
julia> Partition(1, 4)
Domain partitioning with (1, 4, 1) ranks
└── y-partitioning: 4
julia> Partition(x = Fractional(1, 2, 3, 4))
Domain partitioning with (4, 1, 1) ranks
└── x-partitioning: domain fractions: (0.1, 0.2, 0.3, 0.4)
```

`Oceananigans.DistributedComputations.reconstruct_global_grid`

— Method`reconstruct_global_grid(grid::DistributedGrid)`

Return the global grid on `child_architecture(grid)`

## Fields

`Oceananigans.Fields.AbstractField`

— Type`AbstractField{LX, LY, LZ, G, T, N}`

Abstract supertype for fields located at `(LX, LY, LZ)`

and defined on a grid `G`

with eltype `T`

and `N`

dimensions.

Note: we need the parameter `T`

to subtype AbstractArray.

`Oceananigans.Fields.BackgroundField`

— Type`BackgroundField{F, P}`

Temporary container for storing information about `BackgroundFields`

.

`Oceananigans.Fields.BackgroundField`

— Method`BackgroundField(func; parameters=nothing)`

Returns a `BackgroundField`

to be passed to `NonhydrostaticModel`

for use as a background velocity or tracer field.

If `parameters`

is not provided, `func`

must be callable with the signature

`func(x, y, z, t)`

If `parameters`

is provided, `func`

must be callable with the signature

`func(x, y, z, t, parameters)`

`Oceananigans.Fields.Field`

— Method```
Field{LX, LY, LZ}(grid::AbstractGrid,
T::DataType=eltype(grid); kw...) where {LX, LY, LZ}
```

Construct a `Field`

on `grid`

with data type `T`

at the location `(LX, LY, LZ)`

. Each of `(LX, LY, LZ)`

is either `Center`

or `Face`

and determines the field's location in `(x, y, z)`

respectively.

**Keyword arguments**

`data :: OffsetArray`

: An offset array with the fields data. If nothing is provided the field is filled with zeros.`boundary_conditions`

: If nothing is provided, then field is created using the default boundary conditions via`FieldBoundaryConditions`

.`indices`

: Used to prescribe where a reduced field lives on. For example, at which`k`

index does a two-dimensional $x$-$y$ field lives on. Default:`(:, :, :)`

.

**Example**

A field at location `(Face, Face, Center)`

.

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(2, 3, 4), extent=(1, 1, 1));
julia> ω = Field{Face, Face, Center}(grid)
2×3×4 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 8×9×10 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, -2:7) with eltype Float64 with indices -2:5×-2:6×-2:7
└── max=0.0, min=0.0, mean=0.0
```

Now, using `indices`

we can create a two dimensional $x$-$y$ field at location `(Face, Face, Center)`

to compute, e.g., the vertical vorticity $∂v/∂x - ∂u/∂y$ at the fluid's surface $z = 0$, which for `Center`

corresponds to `k = Nz`

.

```
julia> u = XFaceField(grid); v = YFaceField(grid);
julia> ωₛ = Field(∂x(v) - ∂y(u), indices=(:, :, grid.Nz))
2×3×1 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 4:4)
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 8×9×1 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, 4:4) with eltype Float64 with indices -2:5×-2:6×4:4
└── max=0.0, min=0.0, mean=0.0
julia> compute!(ωₛ)
2×3×1 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 4:4)
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 8×9×1 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, 4:4) with eltype Float64 with indices -2:5×-2:6×4:4
└── max=0.0, min=0.0, mean=0.0
```

`Oceananigans.Fields.Reduction`

— Method`Reduction(reduce!, operand; dims)`

Return a `Reduction`

of `operand`

with `reduce!`

, along `dims`

. Note that `Reduction`

expects `reduce!`

to operate in-place.

**Example**

```
using Oceananigans
Nx, Ny, Nz = 3, 3, 3
grid = RectilinearGrid(size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1),
topology=(Periodic, Periodic, Periodic))
c = CenterField(grid)
set!(c, (x, y, z) -> x + y + z)
max_c² = Field(Reduction(maximum!, c^2, dims=3))
compute!(max_c²)
max_c²[1:Nx, 1:Ny]
# output
3×3 Matrix{Float64}:
1.36111 2.25 3.36111
2.25 3.36111 4.69444
3.36111 4.69444 6.25
```

`Oceananigans.Fields.CenterField`

— Function`CenterField(grid, T=eltype(grid); kw...)`

Return a `Field{Center, Center, Center}`

on `grid`

. Additional keyword arguments are passed to the `Field`

constructor.

`Oceananigans.Fields.PressureFields`

— Function`PressureFields(grid, bcs::NamedTuple)`

Return a `NamedTuple`

with pressure fields `pHY′`

and `pNHS`

initialized as `CenterField`

s on `grid`

. Boundary conditions `bcs`

may be specified via a named tuple of `FieldBoundaryCondition`

s.

`Oceananigans.Fields.PressureFields`

— Method`PressureFields(proposed_pressures::NamedTuple{(:pHY′, :pNHS)}, grid, bcs)`

Return a `NamedTuple`

of pressure fields with, overwriting boundary conditions in `proposed_tracer_fields`

with corresponding fields in the `NamedTuple`

`bcs`

.

`Oceananigans.Fields.TendencyFields`

— Method```
TendencyFields(grid, tracer_names;
u = XFaceField(grid),
v = YFaceField(grid),
w = ZFaceField(grid),
kwargs...)
```

Return a `NamedTuple`

with tendencies for all solution fields (velocity fields and tracer fields), initialized on `grid`

. Optional `kwargs`

can be specified to assign data arrays to each tendency field.

`Oceananigans.Fields.TracerFields`

— Method`TracerFields(tracer_names, grid, user_bcs)`

Return a `NamedTuple`

with tracer fields specified by `tracer_names`

initialized as `CenterField`

s on `grid`

. Boundary conditions `user_bcs`

may be specified via a named tuple of `FieldBoundaryCondition`

s.

`Oceananigans.Fields.TracerFields`

— Method`TracerFields(tracer_names, grid; kwargs...)`

Return a `NamedTuple`

with tracer fields specified by `tracer_names`

initialized as `CenterField`

s on `grid`

. Fields may be passed via optional keyword arguments `kwargs`

for each field.

This function is used by `OutputWriters.Checkpointer`

and `TendencyFields`

. ```

`Oceananigans.Fields.TracerFields`

— Method`TracerFields(proposed_tracers::NamedTuple, grid, bcs)`

Return a `NamedTuple`

of tracers, overwriting boundary conditions in `proposed_tracers`

with corresponding fields in the `NamedTuple`

`bcs`

.

`Oceananigans.Fields.TracerFields`

— MethodShortcut constructor for empty tracer fields.

`Oceananigans.Fields.VelocityFields`

— Function`VelocityFields(grid, user_bcs = NamedTuple())`

Return a `NamedTuple`

with fields `u`

, `v`

, `w`

initialized on `grid`

. Boundary conditions `bcs`

may be specified via a named tuple of `FieldBoundaryCondition`

s.

`Oceananigans.Fields.VelocityFields`

— Method`VelocityFields(proposed_velocities::NamedTuple{(:u, :v, :w)}, grid, bcs)`

Return a `NamedTuple`

of velocity fields, overwriting boundary conditions in `proposed_velocities`

with corresponding fields in the `NamedTuple`

`bcs`

.

`Oceananigans.Fields.XFaceField`

— Function`XFaceField(grid, T=eltype(grid); kw...)`

Return a `Field{Face, Center, Center}`

on `grid`

. Additional keyword arguments are passed to the `Field`

constructor.

`Oceananigans.Fields.YFaceField`

— Function`YFaceField(grid, T=eltype(grid); kw...)`

Return a `Field{Center, Face, Center}`

on `grid`

. Additional keyword arguments are passed to the `Field`

constructor.

`Oceananigans.Fields.ZFaceField`

— Function`ZFaceField(grid, T=eltype(grid); kw...)`

Return a `Field{Center, Center, Face}`

on `grid`

. Additional keyword arguments are passed to the `Field`

constructor.

`Oceananigans.Fields.compute!`

— Function`compute!(field)`

Computes `field.data`

from `field.operand`

.

`Oceananigans.Fields.field`

— Method`field(loc, a, grid)`

Build a field from `a`

at `loc`

and on `grid`

.

`Oceananigans.Fields.interior`

— Method`interior(f::Field)`

Returns a view of `f`

that excludes halo points."

`Oceananigans.Fields.interpolate`

— Method`interpolate(at_node, from_field, from_loc, from_grid)`

Interpolate `from_field`

, `at_node`

, on `from_grid`

and at `from_loc`

ation, where `at_node`

is a tuple of coordinates and and `from_loc = (ℓx, ℓy, ℓz)`

.

Note that this is a lower-level `interpolate`

method defined for use in CPU/GPU kernels.

`Oceananigans.Fields.interpolate`

— Method`interpolate(to_node, from_field)`

Interpolate `field`

to the physical point `(x, y, z)`

using trilinear interpolation.

`Oceananigans.Fields.regrid!`

— Method`regrid!(a, b)`

Regrid field `b`

onto the grid of field `a`

.

**Example**

Generate a tracer field on a vertically stretched grid and regrid it on a regular grid.

```
using Oceananigans
Nz, Lz = 2, 1.0
topology = (Flat, Flat, Bounded)
input_grid = RectilinearGrid(size=Nz, z = [0, Lz/3, Lz], topology=topology, halo=1)
input_field = CenterField(input_grid)
input_field[1, 1, 1:Nz] = [2, 3]
output_grid = RectilinearGrid(size=Nz, z=(0, Lz), topology=topology, halo=1)
output_field = CenterField(output_grid)
regrid!(output_field, input_field)
output_field[1, 1, :]
# output
4-element OffsetArray(::Vector{Float64}, 0:3) with eltype Float64 with indices 0:3:
0.0
2.333333333333333
3.0
0.0
```

`Oceananigans.location`

— MethodReturns the location `(LX, LY, LZ)`

of an `AbstractField{LX, LY, LZ}`

.

`Oceananigans.Fields.@compute`

— Macro`@compute(exprs...)`

Call `compute!`

on fields after defining them.

## Forcings

`Oceananigans.Forcings.AdvectiveForcing`

— Method`AdvectiveForcing(u=ZeroField(), v=ZeroField(), w=ZeroField())`

Build a forcing term representing advection by the velocity field `u, v, w`

with an advection `scheme`

.

**Example**

**Using a tracer field to model sinking particles**

```
using Oceananigans
# Physical parameters
gravitational_acceleration = 9.81 # m s⁻²
ocean_density = 1026 # kg m⁻³
mean_particle_density = 2000 # kg m⁻³
mean_particle_radius = 1e-3 # m
ocean_molecular_kinematic_viscosity = 1.05e-6 # m² s⁻¹
# Terminal velocity of a sphere in viscous flow
Δb = gravitational_acceleration * (mean_particle_density - ocean_density) / ocean_density
ν = ocean_molecular_kinematic_viscosity
R = mean_particle_radius
w_Stokes = - 2/9 * Δb / ν * R^2 # m s⁻¹
settling = AdvectiveForcing(w=w_Stokes)
# output
AdvectiveForcing:
├── u: ZeroField{Int64}
├── v: ZeroField{Int64}
└── w: ConstantField(-1.97096)
```

`Oceananigans.Forcings.ContinuousForcing`

— Type`ContinuousForcing{LX, LY, LZ, P, F, D, I, ℑ}`

A callable object that implements a "continuous form" forcing function on a field at the location `LX, LY, LZ`

with optional parameters.

`Oceananigans.Forcings.ContinuousForcing`

— Method`ContinuousForcing(func; parameters=nothing, field_dependencies=())`

Construct a "continuous form" forcing with optional `parameters`

and optional `field_dependencies`

on other fields in a model.

If neither `parameters`

nor `field_dependencies`

are provided, then `func`

must be callable with the signature

`func(X..., t)`

where, on a three-dimensional grid with no `Flat`

directions, `X = (x, y, z)`

is a 3-tuple containing the east-west, north-south, and vertical spatial coordinates, and `t`

is time.

Dimensions with `Flat`

topology are omitted from the coordinate tuple `X`

. For example, on a grid with topology `(Periodic, Periodic, Flat)`

, and with no `parameters`

or `field_dependencies`

, then `func`

must be callable

`func(x, y, t)`

where `x`

and `y`

are the east-west and north-south coordinates, respectively. For another example, on a grid with topology `(Flat, Flat, Bounded)`

(e.g. a single column), and for a forcing with no `parameters`

or `field_dependencies`

, then `func`

must be callable with

`func(z, t)`

where `z`

is the vertical coordinate.

If `field_dependencies`

are provided, the signature of `func`

must include them. For example, if `field_dependencies=(:u, :S)`

(and `parameters`

are *not* provided), and on a three-dimensional grid with no `Flat`

dimensions, then `func`

must be callable with the signature

`func(x, y, z, t, u, S)`

where `u`

is assumed to be the `u`

-velocity component, and `S`

is a tracer. Note that any field which does not have the name `u`

, `v`

, or `w`

is assumed to be a tracer and must be present in `model.tracers`

.

If `parameters`

are provided, then the *last* argument to `func`

must be `parameters`

. For example, if `func`

has no `field_dependencies`

but does depend on `parameters`

, then on a three-dimensional grid it must be callable with the signature

`func(x, y, z, t, parameters)`

With `field_dependencies=(:u, :v, :w, :c)`

and `parameters`

and on a three-dimensional grid, then `func`

must be callable with the signature

`func(x, y, z, t, u, v, w, c, parameters)`

`Oceananigans.Forcings.DiscreteForcing`

— Type`struct DiscreteForcing{P, F}`

Wrapper for "discrete form" forcing functions with optional `parameters`

.

`Oceananigans.Forcings.DiscreteForcing`

— Method`DiscreteForcing(func; parameters=nothing)`

Construct a "discrete form" forcing function with optional parameters. The forcing function is applied at grid point `i, j, k`

.

When `parameters`

are not specified, `func`

must be callable with the signature

`func(i, j, k, grid, clock, model_fields)`

where `grid`

is `model.grid`

, `clock.time`

is the current simulation time and `clock.iteration`

is the current model iteration, and `model_fields`

is a `NamedTuple`

with `u, v, w`

and the fields in `model.tracers`

.

*Note* that the index `end`

does *not* access the final physical grid point of a model field in any direction. The final grid point must be explicitly specified, as in `model_fields.u[i, j, grid.Nz]`

.

When `parameters`

*is* specified, `func`

must be callable with the signature.

`func(i, j, k, grid, clock, model_fields, parameters)`

Above, `parameters`

is, in principle, arbitrary. Note, however, that GPU compilation can place constraints on `typeof(parameters)`

.

`Oceananigans.Forcings.GaussianMask`

— Type`GaussianMask{D}(center, width)`

Callable object that returns a Gaussian masking function centered on `center`

, with `width`

, and varying along direction `D`

, i.e.,

`exp(-(D - center)^2 / (2 * width^2))`

**Examples**

- Create a Gaussian mask centered on
`z=0`

with width`1`

meter.

`julia> mask = GaussianMask{:z}(center=0, width=1)`

`Oceananigans.Forcings.LinearTarget`

— Type`LinearTarget{D}(intercept, gradient)`

Callable object that returns a Linear target function with `intercept`

and `gradient`

, and varying along direction `D`

, i.e.,

`intercept + D * gradient`

**Examples**

Create a linear target function varying in

`z`

, equal to`0`

at`z=0`

and with gradient 10⁻⁶:`julia> target = LinearTarget{:z}(intercept=0, gradient=1e-6)`

`Oceananigans.Forcings.Relaxation`

— Type`struct Relaxation{R, M, T}`

Callable object for restoring fields to a `target`

at some `rate`

and within a `mask`

ed region in `x, y, z`

.

`Oceananigans.Forcings.Relaxation`

— Method`Relaxation(; rate, mask=onefunction, target=zerofunction)`

Returns a `Forcing`

that restores a field to `target(X..., t)`

at the specified `rate`

, in the region `mask(X...)`

.

The functions `onefunction`

and `zerofunction`

always return 1 and 0, respectively. Thus the default `mask`

leaves the whole domain uncovered, and the default `target`

is zero.

**Example**

- Restore a field to zero on a timescale of "3600" (equal to one hour if the time units of the simulation are seconds).

```
using Oceananigans
damping = Relaxation(rate = 1/3600)
# output
Relaxation{Float64, typeof(Oceananigans.Forcings.onefunction), typeof(Oceananigans.Forcings.zerofunction)}
├── rate: 0.0002777777777777778
├── mask: 1
└── target: 0
```

- Restore a field to a linear z-gradient within the bottom 1/4 of a domain on a timescale of "60" (equal to one minute if the time units of the simulation are seconds).

```
dTdz = 0.001 # ⁰C m⁻¹, temperature gradient
T₀ = 20 # ⁰C, surface temperature at z=0
Lz = 100 # m, depth of domain
bottom_sponge_layer = Relaxation(; rate = 1/60,
target = LinearTarget{:z}(intercept=T₀, gradient=dTdz),
mask = GaussianMask{:z}(center=-Lz, width=Lz/4))
# output
Relaxation{Float64, GaussianMask{:z, Float64}, LinearTarget{:z, Float64}}
├── rate: 0.016666666666666666
├── mask: exp(-(z + 100.0)^2 / (2 * 25.0^2))
└── target: 20.0 + 0.001 * z
```

`Oceananigans.Forcings.Forcing`

— Method`Forcing(func; parameters=nothing, field_dependencies=(), discrete_form=false)`

Returns a forcing function added to the tendency of an Oceananigans model field.

If `discrete_form=false`

(the default), and neither `parameters`

nor `field_dependencies`

are provided, then `func`

must be callable with the signature

`func(x, y, z, t)`

where `x, y, z`

are the east-west, north-south, and vertical spatial coordinates, and `t`

is time. Note that this form is also default in the constructor for `NonhydrostaticModel`

, so that `Forcing`

is not needed.

If `discrete_form=false`

(the default), and `field_dependencies`

are provided, the signature of `func`

must include them. For example, if `field_dependencies=(:u, :S)`

(and `parameters`

are *not* provided), then `func`

must be callable with the signature

`func(x, y, z, t, u, S)``

where `u`

is assumed to be the `u`

-velocity component, and `S`

is a tracer. Note that any field which does not have the name `u`

, `v`

, or `w`

is assumed to be a tracer and must be present in `model.tracers`

.

If `discrete_form=false`

(the default) and `parameters`

are provided, then the *last* argument to `func`

must be `parameters`

. For example, if `func`

has no `field_dependencies`

but does depend on `parameters`

, then it must be callable with the signature

`func(x, y, z, t, parameters)`

The object `parameters`

is arbitrary in principle, however GPU compilation can place constraints on `typeof(parameters)`

.

With `field_dependencies=(:u, :v, :w, :c)`

and `parameters`

, then `func`

must be callable with the signature

`func(x, y, z, t, u, v, w, c, parameters)`

If `discrete_form=true`

then `func`

must be callable with the "discrete form"

`func(i, j, k, grid, clock, model_fields)`

where `i, j, k`

is the grid point at which the forcing is applied, `grid`

is `model.grid`

, `clock.time`

is the current simulation time and `clock.iteration`

is the current model iteration, and `model_fields`

is a `NamedTuple`

with `u, v, w`

, the fields in `model.tracers`

, and the fields in `model.diffusivity_fields`

, each of which is an `OffsetArray`

s (or `NamedTuple`

s of `OffsetArray`

s depending on the turbulence closure) of field data.

When `discrete_form=true`

and `parameters`

*is* specified, `func`

must be callable with the signature

`func(i, j, k, grid, clock, model_fields, parameters)`

**Examples**

```
using Oceananigans
# Parameterized forcing
parameterized_func(x, y, z, t, p) = p.μ * exp(z / p.λ) * cos(p.ω * t)
v_forcing = Forcing(parameterized_func, parameters = (μ=42, λ=0.1, ω=π))
# output
ContinuousForcing{NamedTuple{(:μ, :λ, :ω), Tuple{Int64, Float64, Irrational{:π}}}}
├── func: parameterized_func (generic function with 1 method)
├── parameters: (μ = 42, λ = 0.1, ω = π)
└── field dependencies: ()
```

Note that because forcing locations are regularized within the `NonhydrostaticModel`

constructor:

```
grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, forcing=(v=v_forcing,))
model.forcing.v
# output
ContinuousForcing{NamedTuple{(:μ, :λ, :ω), Tuple{Int64, Float64, Irrational{:π}}}} at (Center, Face, Center)
├── func: parameterized_func (generic function with 1 method)
├── parameters: (μ = 42, λ = 0.1, ω = π)
└── field dependencies: ()
```

After passing through the constructor for `NonhydrostaticModel`

, the `v`

-forcing location information is available and set to `Center, Face, Center`

.

```
# Field-dependent forcing
growth_in_sunlight(x, y, z, t, P) = exp(z) * P
plankton_forcing = Forcing(growth_in_sunlight, field_dependencies=:P)
# output
ContinuousForcing{Nothing}
├── func: growth_in_sunlight (generic function with 1 method)
├── parameters: nothing
└── field dependencies: (:P,)
```

```
# Parameterized, field-dependent forcing
tracer_relaxation(x, y, z, t, c, p) = p.μ * exp((z + p.H) / p.λ) * (p.dCdz * z - c)
c_forcing = Forcing(tracer_relaxation,
field_dependencies = :c,
parameters = (μ=1/60, λ=10, H=1000, dCdz=1))
# output
ContinuousForcing{NamedTuple{(:μ, :λ, :H, :dCdz), Tuple{Float64, Int64, Int64, Int64}}}
├── func: tracer_relaxation (generic function with 1 method)
├── parameters: (μ = 0.016666666666666666, λ = 10, H = 1000, dCdz = 1)
└── field dependencies: (:c,)
```

```
# Unparameterized discrete-form forcing function
filtered_relaxation(i, j, k, grid, clock, model_fields) =
@inbounds - (model_fields.c[i-1, j, k] + model_fields.c[i, j, k] + model_fields.c[i+1, j, k]) / 3
filtered_forcing = Forcing(filtered_relaxation, discrete_form=true)
# output
DiscreteForcing{Nothing}
├── func: filtered_relaxation (generic function with 1 method)
└── parameters: nothing
```

```
# Discrete-form forcing function with parameters
masked_damping(i, j, k, grid, clock, model_fields, parameters) =
@inbounds - parameters.μ * exp(grid.zᵃᵃᶜ[k] / parameters.λ) * model_fields.u[i, j, k]
masked_damping_forcing = Forcing(masked_damping, parameters=(μ=42, λ=π), discrete_form=true)
# output
DiscreteForcing{NamedTuple{(:μ, :λ), Tuple{Int64, Irrational{:π}}}}
├── func: masked_damping (generic function with 1 method)
└── parameters: (μ = 42, λ = π)
```

## Grids

`Oceananigans.Grids.AbstractCurvilinearGrid`

— Type`AbstractCurvilinearGrid{FT, TX, TY, TZ}`

Abstract supertype for curvilinear grids with elements of type `FT`

and topology `{TX, TY, TZ}`

.

`Oceananigans.Grids.AbstractGrid`

— Type`AbstractGrid{FT, TX, TY, TZ}`

Abstract supertype for grids with elements of type `FT`

and topology `{TX, TY, TZ}`

.

`Oceananigans.Grids.AbstractHorizontallyCurvilinearGrid`

— Type`AbstractHorizontallyCurvilinearGrid{FT, TX, TY, TZ}`

Abstract supertype for horizontally-curvilinear grids with elements of type `FT`

and topology `{TX, TY, TZ}`

.

`Oceananigans.Grids.AbstractTopology`

— Type`AbstractTopology`

Abstract supertype for grid topologies.

`Oceananigans.Grids.AbstractUnderlyingGrid`

— Type`AbstractUnderlyingGrid{FT, TX, TY, TZ}`

Abstract supertype for "primary" grids (as opposed to grids with immersed boundaries) with elements of type `FT`

and topology `{TX, TY, TZ}`

.

`Oceananigans.Grids.Bounded`

— Type`Bounded`

Grid topology for bounded dimensions, e.g., wall-bounded dimensions.

`Oceananigans.Grids.Center`

— Type`Center`

A type describing the location at the center of a grid cell.

`Oceananigans.Grids.Face`

— Type`Face`

A type describing the location at the face of a grid cell.

`Oceananigans.Grids.Flat`

— Type`Flat`

Grid topology for flat dimensions, generally with one grid point, along which the solution is uniform and does not vary.

`Oceananigans.Grids.FullyConnected`

— Type`FullyConnected`

Grid topology for dimensions that are connected to other models or domains.

`Oceananigans.Grids.LatitudeLongitudeGrid`

— Type```
LatitudeLongitudeGrid([architecture = CPU(), FT = Float64];
size,
longitude,
latitude,
z = nothing,
radius = R_Earth,
topology = nothing,
precompute_metrics = true,
halo = nothing)
```

Creates a `LatitudeLongitudeGrid`

with coordinates `(λ, φ, z)`

denoting longitude, latitude, and vertical coordinate respectively.

**Positional arguments**

`architecture`

: Specifies whether arrays of coordinates and spacings are stored on the CPU or GPU. Default:`CPU()`

.`FT`

: Floating point data type. Default:`Float64`

.

**Keyword arguments**

`size`

(required): A 3-tuple prescribing the number of grid points each direction.`longitude`

(required),`latitude`

(required),`z`

(default:`nothing`

): Each is either a:- 2-tuple that specify the end points of the domain,
- one-dimensional array specifying the cell interface locations, or
- a single-argument function that takes an index and returns cell interface location.

**Note**: the latitude and longitude coordinates extents are expected in degrees.`radius`

: The radius of the sphere the grid lives on. By default is equal to the radius of Earth.`topology`

: Tuple of topologies (`Flat`

,`Bounded`

,`Periodic`

) for each direction. The vertical`topology[3]`

must be`Bounded`

, while the latitude-longitude topologies can be`Bounded`

,`Periodic`

, or`Flat`

. If no topology is provided then, by default, the topology is (`Periodic`

,`Bounded`

,`Bounded`

) if the latitudinal extent is 360 degrees or (`Bounded`

,`Bounded`

,`Bounded`

) otherwise.`precompute_metrics`

: Boolean specifying whether to precompute horizontal spacings and areas. Default:`true`

. When`false`

, horizontal spacings and areas are computed on-the-fly during a simulation.`halo`

: A 3-tuple of integers specifying the size of the halo region of cells surrounding the physical interior. The default is 3 halo cells in every direction.

**Examples**

- A default grid with
`Float64`

type:

```
julia> using Oceananigans
julia> grid = LatitudeLongitudeGrid(size=(36, 34, 25),
longitude = (-180, 180),
latitude = (-85, 85),
z = (-1000, 0))
36×34×25 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [-180.0, 180.0) regularly spaced with Δλ=10.0
├── latitude: Bounded φ ∈ [-85.0, 85.0] regularly spaced with Δφ=5.0
└── z: Bounded z ∈ [-1000.0, 0.0] regularly spaced with Δz=40.0
```

- A bounded spherical sector with cell interfaces stretched hyperbolically near the top:

```
julia> using Oceananigans
julia> σ = 1.1; # stretching factor
julia> Nz = 24; # vertical resolution
julia> Lz = 1000; # depth (m)
julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));
julia> grid = LatitudeLongitudeGrid(size=(36, 34, Nz),
longitude = (-180, 180),
latitude = (-20, 20),
z = hyperbolically_spaced_faces,
topology = (Bounded, Bounded, Bounded))
36×34×24 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── longitude: Bounded λ ∈ [-180.0, 180.0] regularly spaced with Δλ=10.0
├── latitude: Bounded φ ∈ [-20.0, 20.0] regularly spaced with Δφ=1.17647
└── z: Bounded z ∈ [-1000.0, -0.0] variably spaced with min(Δz)=21.3342, max(Δz)=57.2159
```

`Oceananigans.Grids.LeftConnected`

— Type`LeftConnected`

Grid topology for dimensions that are connected to other models or domains only on the left (the other direction is bounded)

`Oceananigans.Grids.Periodic`

— Type`Periodic`

Grid topology for periodic dimensions.

`Oceananigans.Grids.RectilinearGrid`

— Type```
RectilinearGrid([architecture = CPU(), FT = Float64];
size,
x = nothing,
y = nothing,
z = nothing,
halo = nothing,
extent = nothing,
topology = (Periodic, Periodic, Bounded))
```

Create a `RectilinearGrid`

with `size = (Nx, Ny, Nz)`

grid points.

**Positional arguments**

`architecture`

: Specifies whether arrays of coordinates and spacings are stored on the CPU or GPU. Default:`CPU()`

.`FT`

: Floating point data type. Default:`Float64`

.

**Keyword arguments**

`size`

(required): A tuple prescribing the number of grid points in non-`Flat`

directions.`size`

is a 3-tuple for 3D models, a 2-tuple for 2D models, and either a scalar or 1-tuple for 1D models.`topology`

: A 3-tuple`(TX, TY, TZ)`

specifying the topology of the domain.`TX`

,`TY`

, and`TZ`

specify whether the`x`

-,`y`

-, and`z`

directions are`Periodic`

,`Bounded`

, or`Flat`

. The topology`Flat`

indicates that a model does not vary in those directions so that derivatives and interpolation are zero. The default is`topology = (Periodic, Periodic, Bounded)`

.`extent`

: A tuple prescribing the physical extent of the grid in non-`Flat`

directions, e.g.,`(Lx, Ly, Lz)`

. All directions are constructed with regular grid spacing and the domain (in the case that no direction is`Flat`

) is $0 ≤ x ≤ L_x$, $0 ≤ y ≤ L_y$, and $-L_z ≤ z ≤ 0$, which is most appropriate for oceanic applications in which $z = 0$ usually is the ocean's surface.`x`

,`y`

, and`z`

: Each of`x, y, z`

are either (i) 2-tuples that specify the end points of the domain in their respect directions (in which case scalar values may be used in`Flat`

directions), or (ii) arrays or functions of the corresponding indices`i`

,`j`

, or`k`

that specify the locations of cell faces in the`x`

-,`y`

-, or`z`

-direction, respectively. For example, to prescribe the cell faces in`z`

we need to provide a function that takes`k`

as argument and returns the location of the faces for indices`k = 1`

through`k = Nz + 1`

, where`Nz`

is the`size`

of the stretched`z`

dimension.

**Note**: *Either* `extent`

, or *all* of `x`

, `y`

, and `z`

must be specified.

`halo`

: A tuple of integers that specifies the size of the halo region of cells surrounding the physical interior for each non-`Flat`

direction. The default is 3 halo cells in every direction.

The physical extent of the domain can be specified either via `x`

, `y`

, and `z`

keyword arguments indicating the left and right endpoints of each dimensions, e.g., `x = (-π, π)`

or via the `extent`

argument, e.g., `extent = (Lx, Ly, Lz)`

, which specifies the extent of each dimension in which case $0 ≤ x ≤ L_x$, $0 ≤ y ≤ L_y$, and $-L_z ≤ z ≤ 0$.

A grid topology may be specified via a tuple assigning one of `Periodic`

, `Bounded`

, and, `Flat`

to each dimension. By default, a horizontally periodic grid topology `(Periodic, Periodic, Bounded)`

is assumed.

Constants are stored using floating point values of type `FT`

. By default this is `Float64`

. Make sure to specify the desired `FT`

if not using `Float64`

.

**Grid properties**

`(Nx, Ny, Nz) :: Int`

: Number of physical points in the $(x, y, z)$-direction.`(Hx, Hy, Hz) :: Int`

: Number of halo points in the $(x, y, z)$-direction.`(Lx, Ly, Lz) :: FT`

: Physical extent of the grid in the $(x, y, z)$-direction.`(Δxᶜᵃᵃ, Δyᵃᶜᵃ, Δzᵃᵃᶜ)`

: Spacings in the $(x, y, z)$-directions between the cell faces. These are the lengths in $x$, $y$, and $z$ of`Center`

cells and are defined at`Center`

locations.`(Δxᶠᵃᵃ, Δyᵃᶠᵃ, Δzᵃᵃᶠ)`

: Spacings in the $(x, y, z)$-directions between the cell centers. These are the lengths in $x$, $y$, and $z$ of`Face`

cells and are defined at`Face`

locations.`(xᶜᵃᵃ, yᵃᶜᵃ, zᵃᵃᶜ)`

: $(x, y, z)$ coordinates of cell`Center`

s.`(xᶠᵃᵃ, yᵃᶠᵃ, zᵃᵃᶠ)`

: $(x, y, z)$ coordinates of cell`Face`

s.

**Examples**

- A grid with the default
`Float64`

type:

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(32, 32, 32), extent=(1, 2, 3))
32×32×32 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=0.03125
├── Periodic y ∈ [0.0, 2.0) regularly spaced with Δy=0.0625
└── Bounded z ∈ [-3.0, 0.0] regularly spaced with Δz=0.09375
```

- A grid with
`Float32`

type:

```
julia> using Oceananigans
julia> grid = RectilinearGrid(Float32; size=(32, 32, 16), x=(0, 8), y=(-10, 10), z=(-π, π))
32×32×16 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 8.0) regularly spaced with Δx=0.25
├── Periodic y ∈ [-10.0, 10.0) regularly spaced with Δy=0.625
└── Bounded z ∈ [-3.14159, 3.14159] regularly spaced with Δz=0.392699
```

- A two-dimenisional, horizontally-periodic grid:

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(32, 32), extent=(2π, 4π), topology=(Periodic, Periodic, Flat))
32×32×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── Periodic x ∈ [3.60072e-17, 6.28319) regularly spaced with Δx=0.19635
├── Periodic y ∈ [7.20145e-17, 12.5664) regularly spaced with Δy=0.392699
└── Flat z
```

- A one-dimensional "column" grid:

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=256, z=(-128, 0), topology=(Flat, Flat, Bounded))
1×1×256 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded z ∈ [-128.0, 0.0] regularly spaced with Δz=0.5
```

- A horizontally-periodic regular grid with cell interfaces stretched hyperbolically near the top:

```
julia> using Oceananigans
julia> σ = 1.1; # stretching factor
julia> Nz = 24; # vertical resolution
julia> Lz = 32; # depth (m)
julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));
julia> grid = RectilinearGrid(size = (32, 32, Nz),
x = (0, 64), y = (0, 64),
z = hyperbolically_spaced_faces)
32×32×24 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 64.0) regularly spaced with Δx=2.0
├── Periodic y ∈ [0.0, 64.0) regularly spaced with Δy=2.0
└── Bounded z ∈ [-32.0, -0.0] variably spaced with min(Δz)=0.682695, max(Δz)=1.83091
```

- A three-dimensional grid with regular spacing in $x$, cell interfaces at Chebyshev nodes in $y$, and cell interfaces hyperbolically stretched in $z$ near the top:

```
julia> using Oceananigans
julia> Nx, Ny, Nz = 32, 30, 24;
julia> Lx, Ly, Lz = 200, 100, 32; # (m)
julia> chebychev_nodes(j) = - Ly/2 * cos(π * (j - 1) / Ny);
julia> σ = 1.1; # stretching factor
julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));
julia> grid = RectilinearGrid(size = (Nx, Ny, Nz),
topology = (Periodic, Bounded, Bounded),
x = (0, Lx),
y = chebychev_nodes,
z = hyperbolically_spaced_faces)
32×30×24 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 200.0) regularly spaced with Δx=6.25
├── Bounded y ∈ [-50.0, 50.0] variably spaced with min(Δy)=0.273905, max(Δy)=5.22642
└── Bounded z ∈ [-32.0, -0.0] variably spaced with min(Δz)=0.682695, max(Δz)=1.83091
```

`Oceananigans.Grids.RightConnected`

— Type`RightConnected`

Grid topology for dimensions that are connected to other models or domains only on the right (the other direction is bounded)

`Oceananigans.Grids.conformal_cubed_sphere_panel`

— Function```
conformal_cubed_sphere_panel(architecture::AbstractArchitecture = CPU(),
FT::DataType = Float64;
size,
z,
topology = (Bounded, Bounded, Bounded),
ξ = (-1, 1),
η = (-1, 1),
radius = R_Earth,
halo = (1, 1, 1),
rotation = nothing)
```

Create a `OrthogonalSphericalShellGrid`

that represents a section of a sphere after it has been conformally mapped from the face of a cube. The cube's coordinates are `ξ`

and `η`

(which, by default, both take values in the range $[-1, 1]$.

The mapping from the face of the cube to the sphere is done via the CubedSphere.jl package.

**Positional arguments**

`architecture`

: Specifies whether arrays of coordinates and spacings are stored on the CPU or GPU. Default:`CPU()`

.`FT`

: Floating point data type. Default:`Float64`

.

**Keyword arguments**

`size`

(required): A 3-tuple prescribing the number of grid points each direction.`z`

(required): Either a- 2-tuple that specify the end points of the $z$-domain,
- one-dimensional array specifying the cell interface locations, or
- a single-argument function that takes an index and returns cell interface location.

`radius`

: The radius of the sphere the grid lives on. By default this is equal to the radius of Earth.`halo`

: A 3-tuple of integers specifying the size of the halo region of cells surrounding the physical interior. The default is 1 halo cells in every direction.`rotation :: Rotation`

: Rotation of the conformal cubed sphere panel about some axis that passes through the center of the sphere. If`nothing`

is provided (default), then the panel includes the North Pole of the sphere in its center. For example, to construct a grid that includes tha South Pole we can pass either`rotation = RotX(π)`

or`rotation = RotY(π)`

.

**Examples**

- The default conformal cubed sphere panel grid with
`Float64`

type:

```
julia> using Oceananigans, Oceananigans.Grids
julia> grid = conformal_cubed_sphere_panel(size=(36, 34, 25), z=(-1000, 0))
36×34×25 OrthogonalSphericalShellGrid{Float64, Bounded, Bounded, Bounded} on CPU with 1×1×1 halo and with precomputed metrics
├── centered at: North Pole, (λ, φ) = (0.0, 90.0)
├── longitude: Bounded extent 90.0 degrees variably spaced with min(Δλ)=0.616164, max(Δλ)=2.58892
├── latitude: Bounded extent 90.0 degrees variably spaced with min(Δφ)=0.664958, max(Δφ)=2.74119
└── z: Bounded z ∈ [-1000.0, 0.0] regularly spaced with Δz=40.0
```

- The conformal cubed sphere panel that includes the South Pole with
`Float32`

type:

```
julia> using Oceananigans, Oceananigans.Grids, Rotations
julia> grid = conformal_cubed_sphere_panel(Float32, size=(36, 34, 25), z=(-1000, 0), rotation=RotY(π))
36×34×25 OrthogonalSphericalShellGrid{Float32, Bounded, Bounded, Bounded} on CPU with 1×1×1 halo and with precomputed metrics
├── centered at: South Pole, (λ, φ) = (0.0, -90.0)
├── longitude: Bounded extent 90.0 degrees variably spaced with min(Δλ)=0.616167, max(Δλ)=2.58891
├── latitude: Bounded extent 90.0 degrees variably spaced with min(Δφ)=0.664956, max(Δφ)=2.7412
└── z: Bounded z ∈ [-1000.0, 0.0] regularly spaced with Δz=40.0
```

`Oceananigans.Grids.halo_size`

— Method`halo_size(grid)`

Return a 3-tuple with the number of halo cells on either side of the domain in (x, y, z).

`Oceananigans.Grids.minimum_xspacing`

— Method```
minimum_xspacing(grid, ℓx, ℓy, ℓz)
minimum_xspacing(grid) = minimum_xspacing(grid, Center(), Center(), Center())
```

Return the minimum spacing for `grid`

in $x$ direction at location `ℓx, ℓy, ℓz`

.

**Examples**

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));
julia> minimum_xspacing(grid, Center(), Center(), Center())
0.5
```

`Oceananigans.Grids.minimum_yspacing`

— Method```
minimum_yspacing(grid, ℓx, ℓy, ℓz)
minimum_yspacing(grid) = minimum_yspacing(grid, Center(), Center(), Center())
```

Return the minimum spacing for `grid`

in $y$ direction at location `ℓx, ℓy, ℓz`

.

**Examples**

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));
julia> minimum_yspacing(grid, Center(), Center(), Center())
0.25
```

`Oceananigans.Grids.minimum_zspacing`

— Method```
minimum_zspacing(grid, ℓx, ℓy, ℓz)
minimum_zspacing(grid) = minimum_zspacing(grid, Center(), Center(), Center())
```

Return the minimum spacing for `grid`

in $z$ direction at location `ℓx, ℓy, ℓz`

.

**Examples**

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));
julia> minimum_zspacing(grid, Center(), Center(), Center())
0.125
```

`Oceananigans.Grids.new_data`

— Function`new_data(FT, arch, loc, topo, sz, halo_sz, indices)`

Return an `OffsetArray`

of zeros of float type `FT`

on `arch`

itecture, with indices corresponding to a field on a `grid`

of `size(grid)`

and located at `loc`

.

`Oceananigans.Grids.nodes`

— Method```
nodes(grid, (ℓx, ℓy, ℓz); reshape=false, with_halos=false)
nodes(grid, ℓx, ℓy, ℓz; reshape=false, with_halos=false)
```

Return a 3-tuple of views over the interior nodes of the `grid`

's native coordinates at the locations in `loc=(ℓx, ℓy, ℓz)`

in `x, y, z`

.

If `reshape=true`

, the views are reshaped to 3D arrays with non-singleton dimensions 1, 2, 3 for `x, y, z`

, respectively. These reshaped arrays can then be used in broadcast operations with 3D fields or arrays.

For `RectilinearGrid`

s the native coordinates are `x, y, z`

; for curvilinear grids, like `LatitudeLongitudeGrid`

or `OrthogonalSphericalShellGrid`

the native coordinates are `λ, φ, z`

.

`Oceananigans.Grids.offset_data`

— Function`offset_data(underlying_data, grid::AbstractGrid, loc, indices=default_indices(length(loc)))`

Return an `OffsetArray`

that maps to `underlying_data`

in memory, with offset indices appropriate for the `data`

of a field on a `grid`

of `size(grid)`

and located at `loc`

.

`Oceananigans.Grids.topology`

— Method`topology(grid, dim)`

Return the topology of the `grid`

for the `dim`

-th dimension.

`Oceananigans.Grids.topology`

— Method`topology(grid)`

Return a tuple with the topology of the `grid`

for each dimension.

`Oceananigans.Grids.total_size`

— Function`total_size(grid, loc)`

Return the "total" size of a `grid`

at `loc`

. This is a 3-tuple of integers corresponding to the number of grid points along `x, y, z`

.

`Oceananigans.Grids.xnodes`

— Method`xnodes(grid, ℓx, ℓy, ℓz, with_halos=false)`

Return the positions over the interior nodes on `grid`

in the $x$-direction for the location `ℓx`

, `ℓy`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

See `znodes`

for examples.

`Oceananigans.Grids.xspacings`

— Method`xspacings(grid, ℓx, ℓy, ℓz; with_halos=true)`

Return the spacings over the interior nodes on `grid`

in the $x$-direction for the location `ℓx`

, `ℓy`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

```
julia> using Oceananigans
julia> grid = LatitudeLongitudeGrid(size=(8, 15, 10), longitude=(-20, 60), latitude=(-10, 50), z=(-100, 0));
julia> xspacings(grid, Center(), Face(), Center())
16-element view(OffsetArray(::Vector{Float64}, -2:18), 1:16) with eltype Float64:
1.0950562585518518e6
1.1058578920188267e6
1.1112718969963323e6
1.1112718969963323e6
1.1058578920188267e6
1.0950562585518518e6
1.0789196210678827e6
1.0575265956426917e6
1.0309814069457315e6
999413.38046802
962976.3124613502
921847.720658409
876227.979424229
826339.3435524226
772424.8654621692
714747.2110712599
```

`Oceananigans.Grids.ynodes`

— Method`ynodes(grid, ℓx, ℓy, ℓz, with_halos=false)`

Return the positions over the interior nodes on `grid`

in the $y$-direction for the location `ℓx`

, `ℓy`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

See `znodes`

for examples.

`Oceananigans.Grids.yspacings`

— Method`yspacings(grid, ℓx, ℓy, ℓz; with_halos=true)`

Return the spacings over the interior nodes on `grid`

in the $y$-direction for the location `ℓx`

, `ℓy`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

```
julia> using Oceananigans
julia> grid = LatitudeLongitudeGrid(size=(20, 15, 10), longitude=(0, 20), latitude=(-15, 15), z=(-100, 0));
julia> yspacings(grid, Center(), Center(), Center())
222389.85328911748
```

`Oceananigans.Grids.znodes`

— Method`znodes(grid, ℓx, ℓy, ℓz; with_halos=false)`

Return the positions over the interior nodes on `grid`

in the $z$-direction for the location `ℓx`

, `ℓy`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

```
julia> using Oceananigans
julia> horz_periodic_grid = RectilinearGrid(size=(3, 3, 3), extent=(2π, 2π, 1), halo=(1, 1, 1),
topology=(Periodic, Periodic, Bounded));
julia> zC = znodes(horz_periodic_grid, Center())
3-element view(OffsetArray(::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 0:4), 1:3) with eltype Float64:
-0.8333333333333334
-0.5
-0.16666666666666666
julia> zC = znodes(horz_periodic_grid, Center(), Center(), Center())
3-element view(OffsetArray(::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 0:4), 1:3) with eltype Float64:
-0.8333333333333334
-0.5
-0.16666666666666666
julia> zC = znodes(horz_periodic_grid, Center(), Center(), Center(), with_halos=true)
-1.1666666666666667:0.3333333333333333:0.16666666666666666 with indices 0:4
```

`Oceananigans.Grids.zspacings`

— Method`zspacings(grid, ℓx, ℓy, ℓz; with_halos=true)`

Return the spacings over the interior nodes on `grid`

in the $z$-direction for the location `ℓx`

, `ℓy`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

```
julia> using Oceananigans
julia> grid = LatitudeLongitudeGrid(size=(20, 15, 10), longitude=(0, 20), latitude=(-15, 15), z=(-100, 0));
julia> zspacings(grid, Center(), Center(), Center())
10.0
```

`Oceananigans.Grids.λnodes`

— Method`λnodes(grid::AbstractCurvilinearGrid, ℓx, ℓy, ℓz, with_halos=false)`

Return the positions over the interior nodes on a curvilinear `grid`

in the $λ$-direction for the location `ℓλ`

, `ℓφ`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

See `znodes`

for examples.

`Oceananigans.Grids.φnodes`

— Method`φnodes(grid::AbstractCurvilinearGrid, ℓx, ℓy, ℓz, with_halos=false)`

Return the positions over the interior nodes on a curvilinear `grid`

in the $φ$-direction for the location `ℓλ`

, `ℓφ`

, `ℓz`

. For `Bounded`

directions, `Face`

nodes include the boundary points.

See `znodes`

for examples.

## Immersed boundaries

`Oceananigans.ImmersedBoundaries.GridFittedBottom`

— Type`GridFittedBottom(bottom_height, [immersed_condition=CenterImmersedCondition()])`

Return an immersed boundary with an irregular bottom fit to the underlying grid.

`Oceananigans.ImmersedBoundaries.ImmersedBoundaryCondition`

— Method`ImmersedBoundaryCondition(; interfaces...)`

Return an ImmersedBoundaryCondition with conditions on individual cell `interfaces ∈ (west, east, south, north, bottom, top)`

between the fluid and immersed boundary.

`Oceananigans.ImmersedBoundaries.ImmersedBoundaryGrid`

— Method`ImmersedBoundaryGrid(grid, ib::GridFittedBottom)`

Return a grid with `GridFittedBottom`

immersed boundary.

Computes ib.bottom_height and wraps in an array.

## Logger

`Oceananigans.Logger.OceananigansLogger`

— Type`OceananigansLogger(stream::IO=stdout, level=Logging.Info; show_info_source=false)`

Based on Logging.SimpleLogger, it tries to log all messages in the following format:

`[yyyy/mm/dd HH:MM:SS.sss] log_level message [-@-> source_file:line_number]`

where the source of the message between the square brackets is included only if `show_info_source=true`

or if the message is not an info level message.

## Models

`Oceananigans.Models.seawater_density`

— MethodReturn a `KernelFunctionOperation`

to compute the in-situ `seawater_density`

.

`Oceananigans.Models.seawater_density`

— Method`seawater_density(model; temperature, salinity, geopotential_height)`

Return a `KernelFunctionOperation`

that computes the in-situ density of seawater with (gridded) `temperature`

, `salinity`

, and at `geopotential_height`

. To compute the in-situ density, the 55 term polynomial approximation to the equation of state from Roquet et al. (2015) is used. By default the `seawater_density`

extracts the geopotential height from the model to compute the in-situ density. To compute a potential density at some user chosen reference geopotential height, set `geopotential_height`

to a constant for the density computation,

```
geopotential_height = 0 # sea-surface height
σ₀ = seawater_density(model; geopotential_height)
```

**Note:** `seawater_density`

must be passed a `BoussinesqEquationOfState`

to compute the density. See the relevant documentation for how to set `SeawaterBuoyancy`

using a `BoussinesqEquationOfState`

.

**Example**

Compute a density `Field`

using the `KernelFunctionOperation`

returned from `seawater_density`

```
julia> using Oceananigans, SeawaterPolynomials.TEOS10
julia> using Oceananigans.Models: seawater_density
julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded))
1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded z ∈ [-1000.0, 0.0] regularly spaced with Δz=10.0
julia> tracers = (:T, :S)
(:T, :S)
julia> eos = TEOS10EquationOfState()
BoussinesqEquationOfState{Float64}:
├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64}
└── reference_density: 1020.0
julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}
julia> model = NonhydrostaticModel(; grid, buoyancy, tracers)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: Nothing
julia> set!(model, S = 34.7, T = 0.5)
julia> density_operation = seawater_density(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: ρ (generic function with 2 methods)
└── arguments: ("BoussinesqEquationOfState{Float64}", "1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU", "1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU", "KernelFunctionOperation at (Center, Center, Center)")
julia> density_field = Field(density_operation)
1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
└── max=0.0, min=0.0, mean=0.0
julia> compute!(density_field)
1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
└── max=1032.38, min=1027.73, mean=1030.06
```

Values for `temperature`

, `salinity`

and `geopotential_height`

can be passed to `seawater_density`

to override the defaults that are obtained from the `model`

.

### Non-hydrostatic models

`Oceananigans.Models.NonhydrostaticModels.NonhydrostaticModel`

— Method```
NonhydrostaticModel(; grid,
clock = Clock{eltype(grid)}(0, 0, 1),
advection = CenteredSecondOrder(),
buoyancy = nothing,
coriolis = nothing,
stokes_drift = nothing,
forcing::NamedTuple = NamedTuple(),
closure = nothing,
boundary_conditions::NamedTuple = NamedTuple(),
tracers = (),
timestepper = :QuasiAdamsBashforth2,
background_fields::NamedTuple = NamedTuple(),
particles::ParticlesOrNothing = nothing,
```

biogeochemistry::AbstractBGCOrNothing = nothing, velocities = nothing, pressures = nothing, diffusivity*fields = nothing, pressure*solver = nothing, immersed*boundary = nothing, auxiliary*fields = NamedTuple())

Construct a model for a non-hydrostatic, incompressible fluid on `grid`

, using the Boussinesq approximation when `buoyancy != nothing`

. By default, all Bounded directions are rigid and impenetrable.

**Keyword arguments**

`grid`

: (required) The resolution and discrete geometry on which the`model`

is solved. The architecture (CPU/GPU) that the model is solved on is inferred from the architecture of the`grid`

. Note that the grid needs to be regularly spaced in the horizontal dimensions, $x$ and $y$.`advection`

: The scheme that advects velocities and tracers. See`Oceananigans.Advection`

.`buoyancy`

: The buoyancy model. See`Oceananigans.BuoyancyModels`

.`coriolis`

: Parameters for the background rotation rate of the model.`stokes_drift`

: Parameters for Stokes drift fields associated with surface waves. Default:`nothing`

.`forcing`

:`NamedTuple`

of user-defined forcing functions that contribute to solution tendencies.`closure`

: The turbulence closure for`model`

. See`Oceananigans.TurbulenceClosures`

.`boundary_conditions`

:`NamedTuple`

containing field boundary conditions.`tracers`

: A tuple of symbols defining the names of the modeled tracers, or a`NamedTuple`

of preallocated`CenterField`

s.`timestepper`

: A symbol that specifies the time-stepping method. Either`:QuasiAdamsBashforth2`

or`:RungeKutta3`

.`background_fields`

:`NamedTuple`

with background fields (e.g., background flow). Default:`nothing`

.`particles`

: Lagrangian particles to be advected with the flow. Default:`nothing`

.`biogeochemistry`

: Biogeochemical model for`tracers`

.`velocities`

: The model velocities. Default:`nothing`

.`pressures`

: Hydrostatic and non-hydrostatic pressure fields. Default:`nothing`

.`diffusivity_fields`

: Diffusivity fields. Default:`nothing`

.`pressure_solver`

: Pressure solver to be used in the model. If`nothing`

(default), the model constructor chooses the default based on the`grid`

provide.`immersed_boundary`

: The immersed boundary. Default:`nothing`

.`auxiliary_fields`

:`NamedTuple`

of auxiliary fields. Default:`nothing`

### Hydrostatic free-surface models

`Oceananigans.Models.HydrostaticFreeSurfaceModels.ExplicitFreeSurface`

— Type`struct ExplicitFreeSurface{E, T}`

The explicit free surface solver.

`η::Any`

: free surface elevation`gravitational_acceleration::Any`

: gravitational accelerations

`Oceananigans.Models.HydrostaticFreeSurfaceModels.HydrostaticFreeSurfaceModel`

— Method```
HydrostaticFreeSurfaceModel(; grid,
clock = Clock{eltype(grid)}(0, 0, 1),
momentum_advection = CenteredSecondOrder(),
tracer_advection = CenteredSecondOrder(),
buoyancy = SeawaterBuoyancy(eltype(grid)),
coriolis = nothing,
free_surface = ImplicitFreeSurface(gravitational_acceleration = g_Earth),
forcing::NamedTuple = NamedTuple(),
closure = nothing,
boundary_conditions::NamedTuple = NamedTuple(),
tracers = (:T, :S),
particles::ParticlesOrNothing = nothing,
biogeochemistry::AbstractBGCOrNothing = nothing,
velocities = nothing,
pressure = nothing,
diffusivity_fields = nothing,
auxiliary_fields = NamedTuple(),
)
```

Construct a hydrostatic model with a free surface on `grid`

.

**Keyword arguments**

`grid`

: (required) The resolution and discrete geometry on which`model`

is solved. The architecture (CPU/GPU) that the model is solve is inferred from the architecture of the grid.`momentum_advection`

: The scheme that advects velocities. See`Oceananigans.Advection`

.`tracer_advection`

: The scheme that advects tracers. See`Oceananigans.Advection`

.`buoyancy`

: The buoyancy model. See`Oceananigans.BuoyancyModels`

.`coriolis`

: Parameters for the background rotation rate of the model.`forcing`

:`NamedTuple`

of user-defined forcing functions that contribute to solution tendencies.`free_surface`

: The free surface model.`closure`

: The turbulence closure for`model`

. See`Oceananigans.TurbulenceClosures`

.`boundary_conditions`

:`NamedTuple`

containing field boundary conditions.`tracers`

: A tuple of symbols defining the names of the modeled tracers, or a`NamedTuple`

of preallocated`CenterField`

s.`particles`

: Lagrangian particles to be advected with the flow. Default:`nothing`

.`biogeochemistry`

: Biogeochemical model for`tracers`

.`velocities`

: The model velocities. Default:`nothing`

.`pressure`

: Hydrostatic pressure field. Default:`nothing`

.`diffusivity_fields`

: Diffusivity fields. Default:`nothing`

.`auxiliary_fields`

:`NamedTuple`

of auxiliary fields. Default:`nothing`

.

`Oceananigans.Models.HydrostaticFreeSurfaceModels.ImplicitFreeSurface`

— Method`ImplicitFreeSurface(; solver_method=:Default, gravitational_acceleration=g_Earth, solver_settings...)`

Return an implicit free-surface solver. The implicit free-surface equation is

\[\left [ 𝛁_h ⋅ (H 𝛁_h) - \frac{1}{g Δt^2} \right ] η^{n+1} = \frac{𝛁_h ⋅ 𝐐_⋆}{g Δt} - \frac{η^{n}}{g Δt^2} ,\]

where $η^n$ is the free-surface elevation at the $n$-th time step, $H$ is depth, $g$ is the gravitational acceleration, $Δt$ is the time step, $𝛁_h$ is the horizontal gradient operator, and $𝐐_⋆$ is the barotropic volume flux associated with the predictor velocity field $𝐮_⋆$, i.e.,

\[𝐐_⋆ = \int_{-H}^0 𝐮_⋆ \, 𝖽 z ,\]

where

\[𝐮_⋆ = 𝐮^n + \int_{t_n}^{t_{n+1}} 𝐆ᵤ \, 𝖽t .\]

This equation can be solved, in general, using the `PreconditionedConjugateGradientSolver`

but other solvers can be invoked in special cases.

If $H$ is constant, we divide through out to obtain

\[\left ( ∇^2_h - \frac{1}{g H Δt^2} \right ) η^{n+1} = \frac{1}{g H Δt} \left ( 𝛁_h ⋅ 𝐐_⋆ - \frac{η^{n}}{Δt} \right ) .\]

Thus, for constant $H$ and on grids with regular spacing in $x$ and $y$ directions, the free surface can be obtained using the `FFTBasedPoissonSolver`

.

`solver_method`

can be either of:

`:FastFourierTransform`

for`FFTBasedPoissonSolver`

`:HeptadiagonalIterativeSolver`

for`HeptadiagonalIterativeSolver`

`:PreconditionedConjugateGradient`

for`PreconditionedConjugateGradientSolver`

By default, if the grid has regular spacing in the horizontal directions then the `:FastFourierTransform`

is chosen, otherwise the `:HeptadiagonalIterativeSolver`

.

`Oceananigans.Models.HydrostaticFreeSurfaceModels.PrescribedVelocityFields`

— Method```
PrescribedVelocityFields(; u = ZeroField(),
v = ZeroField(),
w = ZeroField(),
parameters = nothing)
```

Builds `PrescribedVelocityFields`

with prescribed functions `u`

, `v`

, and `w`

.

If `isnothing(parameters)`

, then `u, v, w`

are called with the signature

`u(x, y, z, t) = # something interesting`

If `!isnothing(parameters)`

, then `u, v, w`

are called with the signature

`u(x, y, z, t, parameters) = # something parameterized and interesting`

In the constructor for `HydrostaticFreeSurfaceModel`

, the functions `u, v, w`

are wrapped in `FunctionField`

and associated with the model's `grid`

and `clock`

.

`Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitFreeSurface`

— Type`SplitExplicitFreeSurface(; gravitational_acceleration = g_Earth, kwargs...)`

Return a `SplitExplicitFreeSurface`

representing an explicit time discretization of oceanic free surface dynamics with `gravitational_acceleration`

.

**Keyword Arguments**

`substeps`

: The number of substeps that divide the range`(t, t + 2Δt)`

, where`Δt`

is the baroclinic timestep. Note that some averaging functions do not require substepping until`2Δt`

. The number of substeps is reduced automatically to the last index of`averaging_weights`

for which`averaging_weights > 0`

.`cfl`

: If set then the number of`substeps`

are computed based on the advective timescale imposed from the barotropic gravity-wave speed, computed with depth`grid.Lz`

. If`fixed_Δt`

is provided then the number of`substeps`

will adapt to maintain an exact cfl. If not the effective cfl will be always lower than the specified`cfl`

provided that the baroclinic time step`Δt_baroclinic < fixed_Δt`

Either `substeps`

*or* `cfl`

(with `grid`

) need to be prescribed.

`grid`

: Used to compute the corresponding barotropic surface wave speed.`fixed_Δt`

: The maximum baroclinic timestep allowed. If`fixed_Δt`

is a`nothing`

and a cfl is provided, then the number of substeps will be computed on the fly from the baroclinic time step to maintain a constant cfl.`gravitational_acceleration`

: the gravitational acceleration (default:`g_Earth`

)`averaging_kernel`

: function of`τ`

used to average the barotropic transport`U`

and free surface`η`

within the barotropic advancement.`τ`

is the fractional substep going from 0 to 2 with the baroclinic time step`t + Δt`

located at`τ = 1`

. This function should be centered at`τ = 1`

, that is, $∑ (aₘ m /M) = 1$. By default the averaging kernel described by Shchepetkin and McWilliams (2005): https://doi.org/10.1016/j.ocemod.2004.08.002 is chosen.`timestepper`

: Time stepping scheme used for the barotropic advancement. Choose one of:`ForwardBackwardScheme()`

(default):`η = f(U)`

then`U = f(η)`

,`AdamsBashforth3Scheme()`

:`η = f(U, Uᵐ⁻¹, Uᵐ⁻²)`

then`U = f(η, ηᵐ, ηᵐ⁻¹, ηᵐ⁻²)`

.

`Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitFreeSurface`

— Type`struct SplitExplicitFreeSurface`

The split-explicit free surface solver.

`η`

: The instantaneous free surface (`ReducedField`

)`state`

: The entire state for the split-explicit solver (`SplitExplicitState`

)`auxiliary`

: Parameters for timestepping split-explicit solver (`NamedTuple`

)`gravitational_acceleration`

: Gravitational acceleration`settings`

: Settings for the split-explicit scheme (`NamedTuple`

)

### Shallow-water models

`Oceananigans.Models.ShallowWaterModels.ShallowWaterModel`

— Method```
ShallowWaterModel(; grid,
gravitational_acceleration,
clock = Clock{eltype(grid)}(0, 0, 1),
momentum_advection = UpwindBiasedFifthOrder(),
tracer_advection = WENO(),
mass_advection = WENO(),
coriolis = nothing,
forcing::NamedTuple = NamedTuple(),
closure = nothing,
bathymetry = nothing,
tracers = (),
diffusivity_fields = nothing,
boundary_conditions::NamedTuple = NamedTuple(),
timestepper::Symbol = :RungeKutta3,
formulation = ConservativeFormulation())
```

Construct a shallow water model on `grid`

with `gravitational_acceleration`

constant.

**Keyword arguments**

`grid`

: (required) The resolution and discrete geometry on which`model`

is solved. The architecture (CPU/GPU) that the model is solve is inferred from the architecture of the grid.`gravitational_acceleration`

: (required) The gravitational acceleration constant.`clock`

: The`clock`

for the model.`momentum_advection`

: The scheme that advects velocities. See`Oceananigans.Advection`

. Default:`UpwindBiasedFifthOrder()`

.`tracer_advection`

: The scheme that advects tracers. See`Oceananigans.Advection`

. Default:`WENO()`

.`mass_advection`

: The scheme that advects the mass equation. See`Oceananigans.Advection`

. Default:`WENO()`

.`coriolis`

: Parameters for the background rotation rate of the model.`forcing`

:`NamedTuple`

of user-defined forcing functions that contribute to solution tendencies.`closure`

: The turbulence closure for`model`

. See`Oceananigans.TurbulenceClosures`

.`bathymetry`

: The bottom bathymetry.`tracers`

: A tuple of symbols defining the names of the modeled tracers, or a`NamedTuple`

of preallocated`CenterField`

s.`diffusivity_fields`

: Stores diffusivity fields when the closures require a diffusivity to be calculated at each timestep.`boundary_conditions`

:`NamedTuple`

containing field boundary conditions.`timestepper`

: A symbol that specifies the time-stepping method. Either`:QuasiAdamsBashforth2`

or`:RungeKutta3`

(default).`formulation`

: Whether the dynamics are expressed in conservative form (`ConservativeFormulation()`

; default) or in non-conservative form with a vector-invariant formulation for the non-linear terms (`VectorInvariantFormulation()`

).

The `ConservativeFormulation()`

requires `RectilinearGrid`

. Use `VectorInvariantFormulation()`

with `LatitudeLongitudeGrid`

.

`Oceananigans.Models.ShallowWaterModels.ShallowWaterScalarDiffusivity`

— Type```
ShallowWaterScalarDiffusivity([FT::DataType=Float64;]
ν=0, ξ=0, discrete_form=false)
```

Return a scalar diffusivity for the shallow water model.

The diffusivity for the shallow water model is calculated as `h * ν`

so that we get a viscous term in the form $h^{-1} 𝛁 ⋅ (h ν t)$, where $t$ is the 2D stress tensor plus a trace, i.e., $t = 𝛁𝐮 + (𝛁𝐮)^T - ξ I ⋅ (𝛁 ⋅ 𝐮)$.

With the `VectorInvariantFormulation()`

(that evolves $u$ and $v$) we compute $h^{-1} 𝛁(ν h 𝛁 t)$, while with the `ConservativeFormulation()`

(that evolves $u h$ and $v h$) we compute $𝛁 (ν h 𝛁 t)$.

### Lagrangian particle tracking

`Oceananigans.Models.LagrangianParticleTracking.LagrangianParticles`

— Method`LagrangianParticles(particles::StructArray; restitution=1.0, tracked_fields::NamedTuple=NamedTuple(), dynamics=no_dynamics)`

Construct some `LagrangianParticles`

that can be passed to a model. The `particles`

should be a `StructArray`

and can contain custom fields. The coefficient of restitution for particle-wall collisions is specified by `restitution`

.

A number of `tracked_fields`

may be passed in as a `NamedTuple`

of fields. Each particle will track the value of each field. Each tracked field must have a corresponding particle property. So if `T`

is a tracked field, then `T`

must also be a custom particle property.

`dynamics`

is a function of `(lagrangian_particles, model, Δt)`

that is called prior to advecting particles. `parameters`

can be accessed inside the `dynamics`

function.

`Oceananigans.Models.LagrangianParticleTracking.LagrangianParticles`

— Method`LagrangianParticles(; x, y, z, restitution=1.0, dynamics=no_dynamics, parameters=nothing)`

Construct some `LagrangianParticles`

that can be passed to a model. The particles will have initial locations `x`

, `y`

, and `z`

. The coefficient of restitution for particle-wall collisions is specified by `restitution`

.

`dynamics`

is a function of `(lagrangian_particles, model, Δt)`

that is called prior to advecting particles. `parameters`

can be accessed inside the `dynamics`

function.

## MultiRegion

`Oceananigans.MultiRegion.ConformalCubedSphereGrid`

— Type```
ConformalCubedSphereGrid(arch=CPU(), FT=Float64;
panel_size,
z,
horizontal_direction_halo = 1,
z_halo = horizontal_direction_halo,
horizontal_topology = FullyConnected,
z_topology = Bounded,
radius = R_Earth,
partition = CubedSpherePartition(; R = 1),
devices = nothing)
```

Return a `ConformalCubedSphereGrid`

that comprises of six `conformal_cubed_sphere_panel`

grids; we refer to each of these grids as a "panel". Each panel corresponds to a face of the cube.

The keyword arguments prescribe the properties of each of the panels. Only the topology in the vertical direction can be prescribed and that's done via the `z_topology`

keyword argumet (default: `Bounded`

). Topologies in both horizontal directions for a `ConformalCubedSphereGrid`

are *always* `FullyConnected`

.

Halo size in both horizontal dimensions *must* be equal; this is prescribed via the `horizontal_halo :: Integer`

keyword argument. The number of halo points in the $z$-direction is prescribed by the `z_halo :: Integer`

keyword argument.

The connectivity between the `ConformalCubedSphereGrid`

panels is depicted below.

```
+==========+==========+
∥ ↑ ∥ ↑ ∥
∥ 1W ∥ 1S ∥
∥←3N P5 6W→∥←5E P6 2S→∥
∥ 4N ∥ 4E ∥
∥ ↓ ∥ ↓ ∥
+==========+==========+==========+
∥ ↑ ∥ ↑ ∥
∥ 5W ∥ 5S ∥
∥←1N P3 4W→∥←3E P4 6S→∥
∥ 2N ∥ 2E ∥
∥ ↓ ∥ ↓ ∥
+==========+==========+==========+
∥ ↑ ∥ ↑ ∥
∥ 3W ∥ 3S ∥
∥←5N P1 2W→∥←1E P2 4S→∥
∥ 6N ∥ 6E ∥
∥ ↓ ∥ ↓ ∥
+==========+==========+
```

The North Pole of the sphere lies in the center of panel 3 (P3) and the South Pole in the center of panel 6 (P6).

The `partition`

keyword argument prescribes the partitioning in regions within each panel; see `CubedSpherePartition`

. For example, a `CubedSpherePartition(; R=2)`

implies that each of the panels are partitioned into 2 regions in each dimension; this adds up, e.g., to 24 regions for the whole sphere. In the depiction below, the intra-panel `x, y`

indices are depicted in the center of each region and the overall region index is shown at the bottom right of each region.

```
+==========+==========+==========+==========+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ | ∥ | ∥
∥← (1, 2) →|← (2, 2) →∥← (1, 2) →|← (2, 2) →∥
∥ | ∥ | ∥
∥ ↓ 19 | ↓ 20 ∥ ↓ 23 | ↓ 24 ∥
+-------- P 5 --------+-------- P 6 --------+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ | ∥ | ∥
∥← (1, 1) →|← (2, 1) →∥← (1, 1) →|← (2, 1) →∥
∥ | ∥ | ∥
∥ ↓ 17 | ↓ 18 ∥ ↓ 21 | ↓ 22 ∥
+==========+==========+==========+==========+==========+==========+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ | ∥ | ∥
∥← (1, 2) →|← (2, 2) →∥← (1, 2) →|← (2, 2) →∥
∥ | ∥ | ∥
∥ ↓ 11 | ↓ 12 ∥ ↓ 15 | ↓ 16 ∥
+-------- P 3 --------+-------- P 4 --------+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ | ∥ | ∥
∥← (1, 1) →|← (2, 1) →∥← (1, 1) →|← (2, 1) →∥
∥ | ∥ | ∥
∥ ↓ 9 | ↓ 10 ∥ ↓ 13 | ↓ 14 ∥
+==========+==========+==========+==========+==========+==========+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ | ∥ | ∥
∥← (1, 2) →|← (2, 2) →∥← (1, 2) →|← (2, 2) →∥
∥ | ∥ | ∥
∥ ↓ 3 | ↓ 4 ∥ ↓ 7 | ↓ 8 ∥
+-------- P 1 --------+-------- P 2 --------+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ | ∥ | ∥
∥← (1, 1) →|← (2, 1) →∥← (1, 1) →|← (2, 1) →∥
∥ | ∥ | ∥
∥ ↓ 1 | ↓ 2 ∥ ↓ 5 | ↓ 6 ∥
+==========+==========+==========+==========+
```

Below, we show in detail panels 1 and 2 and the connectivity of each panel.

```
+===============+==============+==============+===============+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ 11W | 9W ∥ 9S | 10S ∥
∥←19N (2, 1) 4W→|←3E (2, 2) 7W→∥←4E (2, 1) 8W→|←7E (2, 2) 13S→∥
∥ 1N | 2N ∥ 5N | 6N ∥
∥ ↓ 3 | ↓ 4 ∥ ↓ 7 | ↓ 8 ∥
+------------- P 1 ------------+------------ P 2 -------------+
∥ ↑ | ↑ ∥ ↑ | ↑ ∥
∥ 3S | 4S ∥ 7S | 8S ∥
∥←20N (1, 1) 2W→|←1E (2, 1) 5W→∥←2E (1, 1) 6W→|←5E (2, 1) 14S→∥
∥ 23N | 24N ∥ 24N | 22N ∥
∥ ↓ 1 | ↓ 2 ∥ ↓ 5 | ↓ 6 ∥
+===============+==============+==============+===============+
```

**Example**

```
using Oceananigans
grid = ConformalCubedSphereGrid(panel_size=(12, 12, 1), z=(-1, 0), radius=1)
```

We can find out all connectivities of the regions of our grid. For example, to determine the connectivites on the South boundary of each region we can call

```
using Oceananigans.MultiRegion: CubedSphereRegionalConnectivity, East, West, South, North, getregion
for region in 1:length(grid); println("panel ", region, ": ", getregion(grid.connectivity.connections, 3).south); end
```

`Oceananigans.MultiRegion.ConformalCubedSphereGrid`

— Type```
ConformalCubedSphereGrid(filepath::AbstractString, arch::AbstractArchitecture=CPU(), FT=Float64;
Nz,
z,
panel_halo = (1, 1, 1),
panel_topology = (FullyConnected, FullyConnected, Bounded),
radius = R_Earth,
devices = nothing)
```

Load a `ConformalCubedSphereGrid`

from `filepath`

.

`Oceananigans.MultiRegion.CubedSpherePartition`

— Method`CubedSpherePartition(; R = 1)`

Return a cubed sphere partition with `R`

partitions in each horizontal dimension of each panel of the sphere.

`Oceananigans.MultiRegion.MultiRegionGrid`

— Method```
MultiRegionGrid(global_grid; partition = XPartition(2),
devices = nothing,
validate = true)
```

Split a `global_grid`

into different regions handled by `devices`

.

**Positional Arguments**

`global_grid`

: the grid to be divided into regions.

**Keyword Arguments**

`partition`

: the partitioning required. The implemented partitioning are`XPartition`

(division along the $x$ direction) and`YPartition`

(division along the $y$ direction).`devices`

: the devices to allocate memory on. If`nothing`

is provided (default) then memorey is allocated on the the`CPU`

. For`GPU`

computation it is possible to specify the total number of GPUs or the specific GPUs to allocate memory on. The number of devices does not need to match the number of regions.`validate :: Boolean`

: Whether to validate`devices`

; defautl:`true`

.

**Example**

```
julia> using Oceananigans
julia> grid = RectilinearGrid(size=(10, 12), extent=(1, 1), topology=(Bounded, Bounded, Flat))
10×12×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo
├── Bounded x ∈ [0.0, 1.0] regularly spaced with Δx=0.1
├── Bounded y ∈ [0.0, 1.0] regularly spaced with Δy=0.0833333
└── Flat z
julia> multi_region_grid = MultiRegionGrid(grid, partition = XPartition(5))
┌ Warning: MultiRegion functionalities are experimental: help the development by reporting bugs or non-implemented features!
└ @ Oceananigans.MultiRegion ~/Research/OC.jl/src/MultiRegion/multi_region_grid.jl:53
MultiRegionGrid{Float64, Bounded, Bounded, Flat} partitioned on CPU():
├── grids: 2×12×1 RectilinearGrid{Float64, RightConnected, Bounded, Flat} on CPU with 3×3×0 halo
├── partitioning: Equal partitioning in X with (5 regions)
└── devices: (CPU(), CPU(), CPU(), CPU(), CPU())
```

## Operators

`Oceananigans.Operators.divᶜᶜᶜ`

— Method`divᶜᶜᶜ(i, j, k, grid, u, v, w)`

Calculate the divergence $𝛁·𝐕$ of a vector field $𝐕 = (u, v, w)$,

`1/V * [δxᶜᵃᵃ(Ax * u) + δxᵃᶜᵃ(Ay * v) + δzᵃᵃᶜ(Az * w)]`

which ends up at the cell centers `ccc`

.

`Oceananigans.Operators.ζ₃ᶠᶠᶜ`

— Method`ζ₃ᶠᶠᶜ(i, j, k, grid, u, v)`

The vertical vorticity associated with horizontal velocities $u$ and $v$.

`Oceananigans.Operators.∇²ᶜᶜᶜ`

— Method`∇²ᶜᶜᶜ(i, j, k, grid, c)`

Calculate the Laplacian of $c$ via

`1/V * [δxᶜᵃᵃ(Ax * ∂xᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * ∂yᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * ∂zᵃᵃᶠ(c))]`

which ends up at the location `ccc`

.

## Output readers

`Oceananigans.OutputReaders.FieldDataset`

— Method```
FieldDataset(filepath;
architecture=CPU(), grid=nothing, backend=InMemory(), metadata_paths=["metadata"])
```

Returns a `Dict`

containing a `FieldTimeSeries`

for each field in the JLD2 file located at `filepath`

. Note that model output **must** have been saved with halos.

**Keyword arguments**

`backend`

: Either`InMemory()`

(default) or`OnDisk()`

. The`InMemory`

backend will

load the data fully in memory as a 4D multi-dimensional array while the `OnDisk()`

backend will lazily load field time snapshots when the `FieldTimeSeries`

is indexed linearly.

`metadata_paths`

: A list of JLD2 paths to look for metadata. By default it looks in`file["metadata"]`

.`grid`

: May be specified to override the grid used in the JLD2 file.

`Oceananigans.OutputReaders.FieldTimeSeries`

— Method```
FieldTimeSeries(path, name;
backend = InMemory(),
grid = nothing,
iterations = nothing,
times = nothing)
```

Return a `FieldTimeSeries`

containing a time-series of the field `name`

load from JLD2 output located at `path`

.

**Keyword arguments**

`backend`

:`InMemory()`

to load data into a 4D array,`OnDisk()`

to lazily load data from disk when indexing into`FieldTimeSeries`

.`grid`

: A grid to associate with the data, in the case that the native grid was not serialized properly.`iterations`

: Iterations to load. Defaults to all iterations found in the file.`times`

: Save times to load, as determined through an approximate floating point comparison to recorded save times. Defaults to times associated with`iterations`

. Takes precedence over`iterations`

if`times`

is specified.

`Oceananigans.OutputReaders.FieldTimeSeries`

— Method`FieldTimeSeries{LX, LY, LZ}(grid::AbstractGrid, times; kwargs...) where {LX, LY, LZ} =`

Construct a `FieldTimeSeries`

on `grid`

and at `times`

.

**Keyword arguments**

`indices`

: spatial indices`backend`

: backend,`InMemory(indices=Colon())`

or`OnDisk()`

`path`

: path to data for`backend = OnDisk()`

`name`

: name of field for`backend = OnDisk()`

## Output writers

`Oceananigans.OutputWriters.AveragedTimeInterval`

— Type`mutable struct AveragedTimeInterval <: AbstractSchedule`

Container for parameters that configure and handle time-averaged output.

`Oceananigans.OutputWriters.AveragedTimeInterval`

— Method`AveragedTimeInterval(interval; window=interval, stride=1)`

Returns a `schedule`

that specifies periodic time-averaging of output. The time `window`

specifies the extent of the time-average, which reoccurs every `interval`

.

`output`

is computed and accumulated into the average every `stride`

iterations during the averaging window. For example, `stride=1`

computs output every iteration, whereas `stride=2`

computes output every other iteration. Time-averages with longer `stride`

s are faster to compute, but less accurate.

The time-average of $a$ is a left Riemann sum corresponding to

\[⟨a⟩ = T⁻¹ \int_{tᵢ-T}^{tᵢ} a \mathrm{d} t \, ,\]

where $⟨a⟩$ is the time-average of $a$, $T$ is the time-window for averaging, and the $tᵢ$ are discrete times separated by the time `interval`

. The $tᵢ$ specify both the end of the averaging window and the time at which output is written.

**Example**

```
using Oceananigans.OutputWriters: AveragedTimeInterval
using Oceananigans.Utils: days
schedule = AveragedTimeInterval(4days, window=2days)
# output
AveragedTimeInterval(window=2 days, stride=1, interval=4 days)
```

An `AveragedTimeInterval`

schedule directs an output writer to time-average its outputs before writing them to disk:

```
using Oceananigans
using Oceananigans.OutputWriters: JLD2OutputWriter
using Oceananigans.Utils: minutes
model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)))
simulation = Simulation(model, Δt=10minutes, stop_time=30days)
simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
filename= "averaged_velocity_data.jld2",
schedule = AveragedTimeInterval(4days, window=2days, stride=2))
# output
JLD2OutputWriter scheduled on TimeInterval(4 days):
├── filepath: ./averaged_velocity_data.jld2
├── 3 outputs: (u, v, w) averaged on AveragedTimeInterval(window=2 days, stride=2, interval=4 days)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
```

`Oceananigans.OutputWriters.Checkpointer`

— Method```
Checkpointer(model; schedule,
dir = ".",
prefix = "checkpoint",
overwrite_existing = false,
verbose = false,
cleanup = false,
properties = [:architecture, :grid, :clock, :coriolis,
:buoyancy, :closure, :timestepper, :particles])
```

Construct a `Checkpointer`

that checkpoints the model to a JLD2 file on `schedule.`

The `model.clock.iteration`

is included in the filename to distinguish between multiple checkpoint files.

To restart or "pickup" a model from a checkpoint, specify `pickup = true`

when calling `run!`

, ensuring that the checkpoint file is in directory `dir`

. See `run!`

for more details.

Note that extra model `properties`

can be safely specified, but removing crucial properties such as `:velocities`

will make restoring from the checkpoint impossible.

The checkpointer attempts to serialize as much of the model to disk as possible, but functions or objects containing functions cannot be serialized at this time.

**Keyword arguments**

`schedule`

(required): Schedule that determines when to checkpoint.`dir`

: Directory to save output to. Default:`"."`

(current working directory).`prefix`

: Descriptive filename prefixed to all output files. Default: "checkpoint".`overwrite_existing`

: Remove existing files if their filenames conflict. Default:`false`

.`verbose`

: Log what the output writer is doing with statistics on compute/write times and file sizes. Default:`false`

.`cleanup`

: Previous checkpoint files will be deleted once a new checkpoint file is written. Default:`false`

.`properties`

: List of model properties to checkpoint. This list must contain`[:grid, :architecture, :timestepper, :particles]`

. Default: [:architecture, :grid, :clock, :coriolis, :buoyancy, :closure, :velocities, :tracers, :timestepper, :particles]

`Oceananigans.OutputWriters.JLD2OutputWriter`

— Method```
JLD2OutputWriter(model, outputs; filename, schedule,
dir = ".",
indices = (:, :, :),
with_halos = false,
array_type = Array{Float64},
max_filesize = Inf,
overwrite_existing = false,
init = noinit,
including = [:grid, :coriolis, :buoyancy, :closure],
verbose = false,
part = 1,
jld2_kw = Dict{Symbol, Any}())
```

Construct a `JLD2OutputWriter`

for an Oceananigans `model`

that writes `label, output`

pairs in `outputs`

to a JLD2 file.

The argument `outputs`

may be a `Dict`

or `NamedTuple`

. The keys of `outputs`

are symbols or strings that "name" output data. The values of `outputs`

are either `AbstractField`

s, objects that are called with the signature `output(model)`

, or `WindowedTimeAverage`

s of `AbstractFields`

s, functions, or callable objects.

**Keyword arguments**

**Filenaming**

`filename`

(required): Descriptive filename.`".jld2"`

is appended to`filename`

in the file path if`filename`

does not end in`".jld2"`

.`dir`

: Directory to save output to. Default:`"."`

(current working directory).

**Output frequency and time-averaging**

`schedule`

(required):`AbstractSchedule`

that determines when output is saved.

**Slicing and type conversion prior to output**

`indices`

: Specifies the indices to write to disk with a`Tuple`

of`Colon`

,`UnitRange`

, or`Int`

elements. Indices must be`Colon`

,`Int`

, or contiguous`UnitRange`

. Defaults to`(:, :, :)`

or "all indices". If`!with_halos`

, halo regions are removed from`indices`

. For example,`indices = (:, :, 1)`

will save xy-slices of the bottom-most index.`with_halos`

(Bool): Whether or not to slice halo regions from fields before writing output. Note, that to postprocess saved output (e.g., compute derivatives, etc) information about the boundary conditions is often crucial. In that case you might need to set`with_halos = true`

.`array_type`

: The array type to which output arrays are converted to prior to saving. Default:`Array{Float64}`

.

**File management**

`max_filesize`

: The writer will stop writing to the output file once the file size exceeds`max_filesize`

, and write to a new one with a consistent naming scheme ending in`part1`

,`part2`

, etc. Defaults to`Inf`

.`overwrite_existing`

: Remove existing files if their filenames conflict. Default:`false`

.

**Output file metadata management**

`init`

: A function of the form`init(file, model)`

that runs when a JLD2 output file is initialized. Default:`noinit(args...) = nothing`

.`including`

: List of model properties to save with every file. Default:`[:grid, :coriolis, :buoyancy, :closure]`

**Miscellaneous keywords**

`verbose`

: Log what the output writer is doing with statistics on compute/write times and file sizes. Default:`false`

.`part`

: The starting part number used if`max_filesize`

is finite. Default: 1.`jld2_kw`

: Dict of kwargs to be passed to`jldopen`

when data is written.

**Example**

Write out 3D fields for $u$, $v$, $w$, and a tracer $c$, along with a horizontal average:

```
using Oceananigans
using Oceananigans.Utils: hour, minute
model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)), tracers=(:c,))
simulation = Simulation(model, Δt=12, stop_time=1hour)
function init_save_some_metadata!(file, model)
file["author"] = "Chim Riggles"
file["parameters/coriolis_parameter"] = 1e-4
file["parameters/density"] = 1027
return nothing
end
c_avg = Field(Average(model.tracers.c, dims=(1, 2)))
# Note that model.velocities is NamedTuple
simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
filename = "some_data.jld2",
schedule = TimeInterval(20minute),
init = init_save_some_metadata!)
# output
JLD2OutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: ./some_data.jld2
├── 3 outputs: (u, v, w)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
```

and a time- and horizontal-average of tracer $c$ every 20 minutes of simulation time to a file called `some_averaged_data.jld2`

```
simulation.output_writers[:avg_c] = JLD2OutputWriter(model, (; c=c_avg),
filename = "some_averaged_data.jld2",
schedule = AveragedTimeInterval(20minute, window=5minute))
# output
JLD2OutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: ./some_averaged_data.jld2
├── 1 outputs: c averaged on AveragedTimeInterval(window=5 minutes, stride=1, interval=20 minutes)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
```

`Oceananigans.OutputWriters.NetCDFOutputWriter`

— Method```
NetCDFOutputWriter(model, outputs; filename, schedule
dir = ".",
array_type = Array{Float64},
indices = nothing,
with_halos = false,
global_attributes = Dict(),
output_attributes = Dict(),
dimensions = Dict(),
overwrite_existing = false,
deflatelevel = 0,
verbose = false)
```

Construct a `NetCDFOutputWriter`

that writes `(label, output)`

pairs in `outputs`

(which should be a `Dict`

) to a NetCDF file, where `label`

is a string that labels the output and `output`

is either a `Field`

(e.g. `model.velocities.u`

) or a function `f(model)`

that returns something to be written to disk. Custom output requires the spatial `dimensions`

(a `Dict`

) to be manually specified (see examples).

**Keyword arguments**

**Filenaming**

`filename`

(required): Descriptive filename.`".nc"`

is appended to`filename`

if`filename`

does not end in`".nc"`

.`dir`

: Directory to save output to.

**Output frequency and time-averaging**

`schedule`

(required):`AbstractSchedule`

that determines when output is saved.

**Slicing and type conversion prior to output**

`indices`

: Tuple of indices of the output variables to include. Default is`(:, :, :)`

, which includes the full fields.`with_halos`

: Boolean defining whether or not to include halos in the outputs. Default:`false`

. Note, that to postprocess saved output (e.g., compute derivatives, etc) information about the boundary conditions is often crucial. In that case you might need to set`with_halos = true`

.`array_type`

: The array type to which output arrays are converted to prior to saving. Default:`Array{Float64}`

.`dimensions`

: A`Dict`

of dimension tuples to apply to outputs (required for function outputs).

**File management**

`overwrite_existing`

: If`false`

,`NetCDFOutputWriter`

will be set to append to`filepath`

. If`true`

,`NetCDFOutputWriter`

will overwrite`filepath`

if it exists or create it if not. Default:`false`

. See NCDatasets.jl documentation for more information about its`mode`

option.`deflatelevel`

: Determines the NetCDF compression level of data (integer 0-9; 0 (default) means no compression and 9 means maximum compression). See NCDatasets.jl documentation for more information.

**Miscellaneous keywords**

`global_attributes`

: Dict of model properties to save with every file. Default:`Dict()`

.`output_attributes`

: Dict of attributes to be saved with each field variable (reasonable defaults are provided for velocities, buoyancy, temperature, and salinity; otherwise`output_attributes`

*must*be user-provided).

**Examples**

Saving the $u$ velocity field and temperature fields, the full 3D fields and surface 2D slices to separate NetCDF files:

```
using Oceananigans
grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, tracers=:c)
simulation = Simulation(model, Δt=12, stop_time=3600)
fields = Dict("u" => model.velocities.u, "c" => model.tracers.c)
simulation.output_writers[:field_writer] =
NetCDFOutputWriter(model, fields, filename="fields.nc", schedule=TimeInterval(60))
# output
NetCDFOutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./fields.nc
├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0)
├── 2 outputs: (c, u)
└── array type: Array{Float64}
```

```
simulation.output_writers[:surface_slice_writer] =
NetCDFOutputWriter(model, fields, filename="surface_xy_slice.nc",
schedule=TimeInterval(60), indices=(:, :, grid.Nz))
# output
NetCDFOutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./surface_xy_slice.nc
├── dimensions: zC(1), zF(1), xC(16), yF(16), xF(16), yC(16), time(0)
├── 2 outputs: (c, u)
└── array type: Array{Float64}
```

```
simulation.output_writers[:averaged_profile_writer] =
NetCDFOutputWriter(model, fields,
filename = "averaged_z_profile.nc",
schedule = AveragedTimeInterval(60, window=20),
indices = (1, 1, :))
# output
NetCDFOutputWriter scheduled on TimeInterval(1 minute):
├── filepath: ./averaged_z_profile.nc
├── dimensions: zC(16), zF(17), xC(1), yF(1), xF(1), yC(1), time(0)
├── 2 outputs: (c, u) averaged on AveragedTimeInterval(window=20 seconds, stride=1, interval=1 minute)
└── array type: Array{Float64}
```

`NetCDFOutputWriter`

also accepts output functions that write scalars and arrays to disk, provided that their `dimensions`

are provided:

```
using Oceananigans
Nx, Ny, Nz = 16, 16, 16
grid = RectilinearGrid(size=(Nx, Ny, Nz), extent=(1, 2, 3))
model = NonhydrostaticModel(grid=grid)
simulation = Simulation(model, Δt=1.25, stop_iteration=3)
f(model) = model.clock.time^2; # scalar output
g(model) = model.clock.time .* exp.(znodes(Center, grid)) # vector/profile output
xC, yF = xnodes(grid, Center()), ynodes(grid, Face())
XC = [xC[i] for i in 1:Nx, j in 1:Ny]
YF = [yF[j] for i in 1:Nx, j in 1:Ny]
h(model) = @. model.clock.time * sin(XC) * cos(YF) # xy slice output
outputs = Dict("scalar" => f, "profile" => g, "slice" => h)
dims = Dict("scalar" => (), "profile" => ("zC",), "slice" => ("xC", "yC"))
output_attributes = Dict(
"scalar" => Dict("long_name" => "Some scalar", "units" => "bananas"),
"profile" => Dict("long_name" => "Some vertical profile", "units" => "watermelons"),
"slice" => Dict("long_name" => "Some slice", "units" => "mushrooms")
)
global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7)
simulation.output_writers[:things] =
NetCDFOutputWriter(model, outputs,
schedule=IterationInterval(1), filename="things.nc", dimensions=dims, verbose=true,
global_attributes=global_attributes, output_attributes=output_attributes)
# output
NetCDFOutputWriter scheduled on IterationInterval(1):
├── filepath: ./things.nc
├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0)
├── 3 outputs: (profile, slice, scalar)
└── array type: Array{Float64}
```

`Oceananigans.OutputWriters.WindowedTimeAverage`

— Type`WindowedTimeAverage(operand, model=nothing; schedule)`

Returns an object for computing running averages of `operand`

over `schedule.window`

and recurring on `schedule.interval`

, where `schedule`

is an `AveragedTimeInterval`

. During the collection period, averages are computed every `schedule.stride`

iteration.

`operand`

may be a `Oceananigans.Field`

or a function that returns an array or scalar.

Calling `wta(model)`

for `wta::WindowedTimeAverage`

object returns `wta.result`

.

## Simulations

`Oceananigans.Simulations.Callback`

— Type```
Callback(func, schedule=IterationInterval(1);
parameters=nothing, callsite=TimeStepCallsite())
```

Return `Callback`

that executes `func`

on `schedule`

at the `callsite`

with optional `parameters`

. By default, `schedule = IterationInterval(1)`

and `callsite = TimeStepCallsite()`

.

If `isnothing(parameters)`

, `func(sim::Simulation)`

is called. Otherwise, `func`

is called via `func(sim::Simulation, parameters)`

.

The `callsite`

determines where `Callback`

is executed. The possible values for `callsite`

are

`TimeStepCallsite()`

: after a time-step.`TendencyCallsite()`

: after tendencies are calculated, but before taking a time-step (useful for modifying tendency calculations).`UpdateStateCallsite()`

: within`update_state!`

, after auxiliary variables have been computed (for multi-stage time-steppers,`update_state!`

may be called multiple times per time-step).

`Oceananigans.Simulations.Simulation`

— Method```
Simulation(model; Δt,
verbose = true,
stop_iteration = Inf,
stop_time = Inf,
wall_time_limit = Inf)
```

Construct a `Simulation`

for a `model`

with time step `Δt`

.

**Keyword arguments**

`Δt`

: Required keyword argument specifying the simulation time step. Can be a`Number`

for constant time steps or a`TimeStepWizard`

for adaptive time-stepping.`stop_iteration`

: Stop the simulation after this many iterations.`stop_time`

: Stop the simulation once this much model clock time has passed.`wall_time_limit`

: Stop the simulation if it's been running for longer than this many seconds of wall clock time.

`Oceananigans.Simulations.TimeStepWizard`

— Type```
TimeStepWizard([FT=Float64;]
cfl = 0.2,
diffusive_cfl = Inf,
max_change = 1.1,
min_change = 0.5,
max_Δt = Inf,
min_Δt = 0.0,
cell_advection_timescale = cell_advection_timescale,
cell_diffusion_timescale = infinite_diffusion_timescale)
```

Callback function that adjusts the simulation time step to meet specified target values for advective and diffusive Courant-Friedrichs-Lewy (CFL) numbers (`cfl`

and `diffusive_cfl`

), subject to the limits

`max(min_Δt, min_change * previous_Δt) ≤ new_Δt ≤ min(max_Δt, max_change * previous_Δt)`

where `new_Δt`

is the new time step calculated by the `TimeStepWizard`

.

For more information on the CFL number, see its wikipedia entry.

**Example**

To use `TimeStepWizard`

, insert it into a `Callback`

and then add the `Callback`

to a `Simulation`

:

```
julia> simulation = Simulation(model, Δt=0.9, stop_iteration=100)
julia> wizard = TimeStepWizard(cfl=0.2)
julia> simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))
```

Then when `run!(simulation)`

is invoked, the time-step `simulation.Δt`

will be updated every 4 iterations.

(Note that the name `:wizard`

is unimportant.)

`Oceananigans.Models.iteration`

— Method`iteration(sim::Simulation)`

Return the current simulation iteration.

`Oceananigans.Simulations.add_callback!`

— Method```
add_callback!(simulation, callback::Callback; name = GenericName(), callback_kw...)
add_callback!(simulation, func, schedule=IterationInterval(1); name = GenericName(), callback_kw...)
```

Add `Callback(func, schedule)`

to `simulation.callbacks`

under `name`

. The default `GenericName()`

generates a name of the form `:callbackN`

, where `N`

is big enough for the name to be unique.

If `name`

is supplied, it may be modified if `simulation.callbacks[name]`

already exists.

`callback_kw`

are passed to the constructor for `Callback`

.

The `callback`

(which contains a schedule) can also be supplied directly.

`Oceananigans.Simulations.conjure_time_step_wizard!`

— Function`conjure_time_step_wizard!(simulation, schedule=IterationInterval(5), wizard_kw...)`

Add a `TimeStepWizard`

built with `wizard_kw`

as a `Callback`

to `simulation`

, called on `schedule`

which is `IterationInterval(5)`

by default.

`Oceananigans.Simulations.run!`

— Method`run!(simulation; pickup=false)`

Run a `simulation`

until one of `simulation.stop_criteria`

evaluates `true`

. The simulation will then stop.

**Picking simulations up from a checkpoint**

Simulations are "picked up" from a checkpoint if `pickup`

is either `true`

, a `String`

, or an `Integer`

greater than 0.

Picking up a simulation sets field and tendency data to the specified checkpoint, leaving all other model properties unchanged.

Possible values for `pickup`

are:

`pickup=true`

picks a simulation up from the latest checkpoint associated with the`Checkpointer`

in`simulation.output_writers`

.`pickup=iteration::Int`

picks a simulation up from the checkpointed file associated with`iteration`

and the`Checkpointer`

in`simulation.output_writers`

.`pickup=filepath::String`

picks a simulation up from checkpointer data in`filepath`

.

Note that `pickup=true`

and `pickup=iteration`

fails if `simulation.output_writers`

contains more than one checkpointer.

## Solvers

`Oceananigans.Solvers.BatchedTridiagonalSolver`

— Type`struct BatchedTridiagonalSolver{A, B, C, T, G, P}`

A batched solver for large numbers of triadiagonal systems.

`Oceananigans.Solvers.BatchedTridiagonalSolver`

— Method```
BatchedTridiagonalSolver(grid;
lower_diagonal,
diagonal,
upper_diagonal,
scratch = arch_array(architecture(grid), zeros(eltype(grid), size(grid)...)),
tridiagonal_direction = ZDirection()
parameters = nothing)
```

Construct a solver for batched tridiagonal systems on `grid`

of the form

```
bⁱʲ¹ ϕⁱʲ¹ + cⁱʲ¹ ϕⁱʲ² = fⁱʲ¹,
aⁱʲᵏ⁻¹ ϕⁱʲᵏ⁻¹ + bⁱʲᵏ ϕⁱʲᵏ + cⁱʲᵏ ϕⁱʲᵏ⁺¹ = fⁱʲᵏ, k = 2, ..., N-1
aⁱʲᴺ⁻¹ ϕⁱʲᴺ⁻¹ + bⁱʲᴺ ϕⁱʲᴺ = fⁱʲᴺ,
```

or in matrix form

```
⎡ bⁱʲ¹ cⁱʲ¹ 0 ⋯ 0 ⎤ ⎡ ϕⁱʲ¹ ⎤ ⎡ fⁱʲ¹ ⎤
⎢ aⁱʲ¹ bⁱʲ² cⁱʲ² 0 ⋯ ⋮ ⎥ ⎢ ϕⁱʲ² ⎥ ⎢ fⁱʲ² ⎥
⎢ 0 ⋱ ⋱ ⋱ ⎥ ⎢ . ⎥ ⎢ . ⎥
⎢ ⋮ 0 ⎥ ⎢ ϕⁱʲᵏ ⎥ ⎢ fⁱʲᵏ ⎥
⎢ ⋮ aⁱʲᴺ⁻² bⁱʲᴺ⁻¹ cⁱʲᴺ⁻¹ ⎥ ⎢ ⎥ ⎢ . ⎥
⎣ 0 ⋯ 0 aⁱʲᴺ⁻¹ bⁱʲᴺ ⎦ ⎣ ϕⁱʲᴺ ⎦ ⎣ fⁱʲᴺ ⎦
```

where `a`

is the `lower_diagonal`

, `b`

is the `diagonal`

, and `c`

is the `upper_diagonal`

.

Note the convention used here for indexing the upper and lower diagonals; this can be different from other implementations where, e.g., `aⁱʲ²`

may appear at the second row, instead of `aⁱʲ¹`

as above.

`ϕ`

is the solution and `f`

is the right hand side source term passed to `solve!(ϕ, tridiagonal_solver, f)`

.

`a`

, `b`

, `c`

, and `f`

can be specified in three ways:

A 1D array means, e.g., that

`aⁱʲᵏ = a[k]`

.A 3D array means, e.g., that

`aⁱʲᵏ = a[i, j, k]`

.

Other coefficient types can be implemented by extending `get_coefficient`

.

`Oceananigans.Solvers.FFTBasedPoissonSolver`

— Type`FFTBasedPoissonSolver(grid, planner_flag=FFTW.PATIENT)`

Return an `FFTBasedPoissonSolver`

that solves the "generalized" Poisson equation,

\[(∇² + m) ϕ = b,\]

where $m$ is a number, using a eigenfunction expansion of the discrete Poisson operator on a staggered grid and for periodic or Neumann boundary conditions.

In-place transforms are applied to $b$, which means $b$ must have complex-valued elements (typically the same type as `solver.storage`

).

See `solve!`

for more information about the FFT-based Poisson solver algorithm.

`Oceananigans.Solvers.HeptadiagonalIterativeSolver`

— Method```
HeptadiagonalIterativeSolver(coeffs;
grid,
iterative_solver = cg!,
maximum_iterations = prod(size(grid)),
tolerance = 1e-13,
reduced_dim = (false, false, false),
placeholder_timestep = -1.0,
preconditioner_method = :Default,
preconditioner_settings = nothing,
template = arch_array(architecture(grid), zeros(prod(size(grid)))),
verbose = false)
```

Return a `HeptadiagonalIterativeSolver`

to solve the problem `A * x = b`

, provided that `A`

is a symmetric matrix.

The solver relies on a sparse version of the matrix `A`

that is stored in `matrix_constructors`

.

In particular, given coefficients `Ax`

, `Ay`

, `Az`

, `C`

, `D`

, the solved problem is

```
Axᵢ₊₁ ηᵢ₊₁ + Axᵢ ηᵢ₋₁ + Ayⱼ₊₁ ηⱼ₊₁ + Ayⱼ ηⱼ₋₁ + Azₖ₊₁ ηₖ₊₁ + Azₖ ηₖ₋₁
- 2 ( Axᵢ₊₁ + Axᵢ + Ayⱼ₊₁ + Ayⱼ + Azₖ₊₁ + Azₖ ) ηᵢⱼₖ
+ ( Cᵢⱼₖ + Dᵢⱼₖ/Δt^2 ) ηᵢⱼₖ = b
```

To have the equation solved at location `{Center, Center, Center}`

, the coefficients must be specified at:

`Ax`

->`{Face, Center, Center}`

`Ay`

->`{Center, Face, Center}`

`Az`

->`{Center, Center, Face}`

`C`

->`{Center, Center, Center}`

`D`

->`{Center, Center, Center}`

`solver.matrix`

is precomputed with a placeholder timestep value of `placeholder_timestep = -1.0`

.

The sparse matrix `A`

can be constructed with:

`SparseMatrixCSC(constructors...)`

for CPU`CuSparseMatrixCSC(constructors...)`

for GPU

The matrix constructors are calculated based on the pentadiagonal coeffients passed as an input to `matrix_from_coefficients`

function.

To allow for variable time step, the diagonal term `- Az / (g * Δt²)`

is only added later on and it is updated only when the previous time step changes (`previous_Δt != Δt`

).

Preconditioning is done through the various methods implemented in `Solvers/sparse_preconditioners.jl`

.

The `iterative_solver`

used can is to be chosen from the IterativeSolvers.jl package. The default solver is a Conjugate Gradient (`cg`

):

`solver = HeptadiagonalIterativeSolver((Ax, Ay, Az, C, D); grid)`

`Oceananigans.Solvers.PreconditionedConjugateGradientSolver`

— Method```
PreconditionedConjugateGradientSolver(linear_operation;
template_field,
maxiter = size(template_field.grid),
reltol = sqrt(eps(eltype(template_field.grid))),
abstol = 0,
preconditioner = nothing)
```

Returns a `PreconditionedConjugateGradientSolver`

that solves the linear equation $A x = b$ using a iterative conjugate gradient method with optional preconditioning.

The solver is used by calling

`solve!(x, solver::PreconditionedConjugateGradientOperator, b, args...)`

for `solver`

, right-hand side `b`

, solution `x`

, and optional arguments `args...`

.

**Arguments**

`linear_operation`

: Function with signature`linear_operation!(p, y, args...)`

that calculates`A * y`

and stores the result in`p`

for a "candidate solution"`y`

.`args...`

are optional positional arguments passed from`solve!(x, solver, b, args...)`

.`template_field`

: Dummy field that is the same type and size as`x`

and`b`

, which is used to infer the`architecture`

,`grid`

, and to create work arrays that are used internally by the solver.`maxiter`

: Maximum number of iterations the solver may perform before exiting.`reltol, abstol`

: Relative and absolute tolerance for convergence of the algorithm. The iteration stops when`norm(A * x - b) < tolerance`

.`preconditioner`

: Object for which`precondition!(z, preconditioner, r, args...)`

computes`z = P * r`

, where`r`

is the residual. Typically`P`

is approximately`A⁻¹`

.

See `solve!`

for more information about the preconditioned conjugate-gradient algorithm.

`Oceananigans.Solvers.solve!`

— Function`solve!(ϕ, solver::FFTBasedPoissonSolver, b, m=0)`

Solves the "generalized" Poisson equation,

\[(∇² + m) ϕ = b,\]

where $m$ is a number, using a eigenfunction expansion of the discrete Poisson operator on a staggered grid and for periodic or Neumann boundary conditions.

In-place transforms are applied to $b$, which means $b$ must have complex-valued elements (typically the same type as `solver.storage`

).

Equation $(∇² + m) ϕ = b$ is sometimes referred to as the "screened Poisson" equation when $m < 0$, or the Helmholtz equation when $m > 0$.

`Oceananigans.Solvers.solve!`

— Method`solve!(ϕ, solver::BatchedTridiagonalSolver, rhs, args...)`

Solve the batched tridiagonal system of linear equations with right hand side `rhs`

and lower diagonal, diagonal, and upper diagonal coefficients described by the `BatchedTridiagonalSolver`

`solver`

. `BatchedTridiagonalSolver`

uses a modified TriDiagonal Matrix Algorithm (TDMA).

The result is stored in `ϕ`

which must have size `(grid.Nx, grid.Ny, grid.Nz)`

.

Implementation follows Press William *et al.* (1992); §2.4. Note that a slightly different notation from Press et al. is used for indexing the off-diagonal elements; see `BatchedTridiagonalSolver`

.

**Reference**

Press William, H., Teukolsky Saul, A., Vetterling William, T., & Flannery Brian, P. (1992). Numerical recipes: the art of scientific computing. Cambridge University Press

`Oceananigans.Solvers.solve!`

— Method`solve!(x, solver::PreconditionedConjugateGradientSolver, b, args...)`

Solve `A * x = b`

using an iterative conjugate-gradient method, where `A * x`

is determined by `solver.linear_operation`

See figure 2.5 in

The Preconditioned Conjugate Gradient Method in "Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods" Barrett et. al, 2nd Edition.

Given:

- Linear Preconditioner operator
`M!(solution, x, other_args...)`

that computes`M * x = solution`

- A matrix operator
`A`

as a function`A()`

; - A dot product function
`norm()`

; - A right-hand side
`b`

; - An initial guess
`x`

; and - Local vectors:
`z`

,`r`

,`p`

,`q`

This function executes the psuedocode algorithm

```
β = 0
r = b - A(x)
iteration = 0
Loop:
if iteration > maxiter
break
end
ρ = r ⋅ z
z = M(r)
β = ρⁱ⁻¹ / ρ
p = z + β * p
q = A(p)
α = ρ / (p ⋅ q)
x = x + α * p
r = r - α * q
if |r| < tolerance
break
end
iteration += 1
ρⁱ⁻¹ = ρ
```

## Stokes drift

`Oceananigans.StokesDrifts.StokesDrift`

— Method```
StokesDrift(; ∂z_uˢ=zerofunction, ∂y_uˢ=zerofunction, ∂t_uˢ=zerofunction,
∂z_vˢ=zerofunction, ∂x_vˢ=zerofunction, ∂t_vˢ=zerofunction,
∂x_wˢ=zerofunction, ∂y_wˢ=zerofunction, ∂t_wˢ=zerofunction, parameters=nothing)
```

Construct a set of functions of space and time for a Stokes drift velocity field corresponding to a surface gravity wave field with an envelope that (potentially) varies in the horizontal directions.

To resolve the evolution of the Lagrangian-mean momentum, we require all the components of the "psuedovorticity",

\[𝛁 × 𝐯ˢ = \hat{\boldsymbol{x}} (∂_y wˢ - ∂_z vˢ) + \hat{\boldsymbol{y}} (∂_z uˢ - ∂_x wˢ) + \hat{\boldsymbol{z}} (∂_x vˢ - ∂_y uˢ)\]

as well as the time-derivatives of $uˢ$, $vˢ$, and $wˢ$.

Note that each function (e.g., `∂z_uˢ`

) is generally a function of depth, horizontal coordinates, and time.Thus, the correct function signature depends on the grid, since `Flat`

horizontal directions are omitted.

For example, on a grid with `topology = (Periodic, Flat, Bounded)`

(and `parameters=nothing`

), then, e.g., `∂z_uˢ`

is callable via `∂z_uˢ(x, z, t)`

. When `!isnothing(parameters)`

, then `∂z_uˢ`

is callable via `∂z_uˢ(x, z, t, parameters)`

. Similarly, on a grid with `topology = (Periodic, Periodic, Bounded)`

and `parameters=nothing`

, `∂z_uˢ`

is called via `∂z_uˢ(x, y, z, t)`

.

**Example**

A wavepacket moving with the group velocity in the $x$-direction. We write the Stokes drift as:

\[uˢ(x, y, z, t) = A(x - cᵍ \, t, y) ûˢ(z)\]

with $A(ξ, η) = \exp{[-(ξ^2 + η^2) / 2δ^2]}$. We also assume $vˢ = 0$. If $𝐯ˢ$ represents the solenoidal component of the Stokes drift, then in this system from incompressibility requirement we have that $∂_z wˢ = - ∂_x uˢ = - (∂_ξ A) ûˢ$ and therefore, under the assumption that $wˢ$ tends to zero at large depths, we get $wˢ = - (∂_ξ A / 2k) ûˢ$.

```
using Oceananigans
using Oceananigans.Units
g = 9.81 # gravitational acceleration
ϵ = 0.1
λ = 100meters # horizontal wavelength
const k = 2π / λ # horizontal wavenumber
c = sqrt(g / k) # phase speed
const δ = 400kilometers # wavepacket spread
const cᵍ = c / 2 # group speed
const Uˢ = ϵ^2 * c
@inline A(ξ, η) = exp(- (ξ^2 + η^2) / 2δ^2)
@inline ∂ξ_A(ξ, η) = - ξ / δ^2 * A(ξ, η)
@inline ∂η_A(ξ, η) = - η / δ^2 * A(ξ, η)
@inline ∂η_∂ξ_A(ξ, η) = η * ξ / δ^4 * A(ξ, η)
@inline ∂²ξ_A(ξ, η) = (ξ^2 / δ^2 - 1) * A(ξ, η) / δ^2
@inline ûˢ(z) = Uˢ * exp(2k * z)
@inline uˢ(x, y, z, t) = A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂z_uˢ(x, y, z, t) = 2k * A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂y_uˢ(x, y, z, t) = ∂η_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂t_uˢ(x, y, z, t) = - cᵍ * ∂ξ_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂x_wˢ(x, y, z, t) = - 1 / 2k * ∂²ξ_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂y_wˢ(x, y, z, t) = - 1 / 2k * ∂η_∂ξ_A(x - cᵍ * t, y) * ûˢ(z)
@inline ∂t_wˢ(x, y, z, t) = + cᵍ / 2k * ∂²ξ_A(x - cᵍ * t, y) * ûˢ(z)
stokes_drift = StokesDrift(; ∂z_uˢ, ∂t_uˢ, ∂y_uˢ, ∂t_wˢ, ∂x_wˢ, ∂y_wˢ)
# output
StokesDrift{Nothing}:
├── ∂x_vˢ: zerofunction
├── ∂x_wˢ: ∂x_wˢ
├── ∂y_uˢ: ∂y_uˢ
├── ∂y_wˢ: ∂y_wˢ
├── ∂z_uˢ: ∂z_uˢ
├── ∂z_vˢ: zerofunction
├── ∂t_uˢ: ∂t_uˢ
├── ∂t_vˢ: zerofunction
└── ∂t_wˢ: ∂t_wˢ
```

`Oceananigans.StokesDrifts.UniformStokesDrift`

— Method`UniformStokesDrift(; ∂z_uˢ=zerofunction, ∂z_vˢ=zerofunction, ∂t_uˢ=zerofunction, ∂t_vˢ=zerofunction, parameters=nothing)`

Construct a set of functions for a Stokes drift velocity field corresponding to a horizontally-uniform surface gravity wave field, with optional `parameters`

.

If `parameters=nothing`

, then the functions `∂z_uˢ`

, `∂z_vˢ`

, `∂t_uˢ`

, `∂t_vˢ`

must be callable with signature `(z, t)`

. If `!isnothing(parameters)`

, then functions must be callable with the signature `(z, t, parameters)`

.

To resolve the evolution of the Lagrangian-mean momentum, we require vertical-derivatives and time-derivatives of the horizontal components of the Stokes drift, `uˢ`

and `vˢ`

.

**Examples**

Exponentially decaying Stokes drift corresponding to a surface Stokes drift of `uˢ(z=0) = 0.005`

and decay scale `h = 20`

:

```
using Oceananigans
@inline uniform_stokes_shear(z, t) = 0.005 * exp(z / 20)
stokes_drift = UniformStokesDrift(∂z_uˢ=uniform_stokes_shear)
# output
UniformStokesDrift{Nothing}:
├── ∂z_uˢ: uniform_stokes_shear
├── ∂z_vˢ: zerofunction
├── ∂t_uˢ: zerofunction
└── ∂t_vˢ: zerofunction
```

Exponentially-decaying Stokes drift corresponding to a surface Stokes drift of `uˢ = 0.005`

and decay scale `h = 20`

, using parameters:

```
using Oceananigans
@inline uniform_stokes_shear(z, t, p) = p.uˢ * exp(z / p.h)
stokes_drift_parameters = (uˢ = 0.005, h = 20)
stokes_drift = UniformStokesDrift(∂z_uˢ=uniform_stokes_shear, parameters=stokes_drift_parameters)
# output
UniformStokesDrift with parameters (uˢ=0.005, h=20):
├── ∂z_uˢ: uniform_stokes_shear
├── ∂z_vˢ: zerofunction
├── ∂t_uˢ: zerofunction
└── ∂t_vˢ: zerofunction
```

## Time steppers

`Oceananigans.TimeSteppers.Clock`

— Type`mutable struct Clock{T<:Number}`

Keeps track of the current `time`

, `iteration`

number, and time-stepping `stage`

. The `stage`

is updated only for multi-stage time-stepping methods. The `time::T`

is either a number or a `DateTime`

object.

`Oceananigans.TimeSteppers.Clock`

— Method`Clock(; time, iteration=0, stage=1)`

Returns a `Clock`

object. By default, `Clock`

is initialized to the zeroth `iteration`

and first time step `stage`

.

`Oceananigans.TimeSteppers.QuasiAdamsBashforth2TimeStepper`

— Method```
QuasiAdamsBashforth2TimeStepper(grid, tracers,
χ = 0.1;
implicit_solver = nothing,
Gⁿ = TendencyFields(grid, tracers),
G⁻ = TendencyFields(grid, tracers))
```

Return a 2nd-order quasi Adams-Bashforth (AB2) time stepper (`QuasiAdamsBashforth2TimeStepper`

) on `grid`

, with `tracers`

, and AB2 parameter `χ`

. The tendency fields `Gⁿ`

and `G⁻`

can be specified via optional `kwargs`

.

The 2nd-order quasi Adams-Bashforth timestepper steps forward the state `Uⁿ`

by `Δt`

via

`Uⁿ⁺¹ = Uⁿ + Δt * [(3/2 + χ) * Gⁿ - (1/2 + χ) * Gⁿ⁻¹]`

where `Uⁿ`

is the state at the $n$-th timestep, `Gⁿ`

is the tendency at the $n$-th timestep, and `Gⁿ⁻¹`

is the tendency at the previous timestep (`G⁻`

).

For the first timestep, since there are no saved tendencies from the previous timestep, the `QuasiAdamsBashforth2TimeStepper`

performs an Euler timestep:

`Uⁿ⁺¹ = Uⁿ + Δt * Gⁿ`

`Oceananigans.TimeSteppers.RungeKutta3TimeStepper`

— Type`RungeKutta3TimeStepper{FT, TG} <: AbstractTimeStepper`

Holds parameters and tendency fields for a low storage, third-order Runge-Kutta-Wray time-stepping scheme described by Le and Moin (1991).

`Oceananigans.TimeSteppers.RungeKutta3TimeStepper`

— Method```
RungeKutta3TimeStepper(grid, tracers;
implicit_solver = nothing,
Gⁿ = TendencyFields(grid, tracers),
G⁻ = TendencyFields(grid, tracers))
```

Return a 3rd-order Runge0Kutta timestepper (`RungeKutta3TimeStepper`

) on `grid`

and with `tracers`

. The tendency fields `Gⁿ`

and `G⁻`

can be specified via optional `kwargs`

.

The scheme described by Le and Moin (1991). In a nutshel, the 3rd-order Runge Kutta timestepper steps forward the state `Uⁿ`

by `Δt`

via 3 substeps. A pressure correction step is applied after at each substep.

The state `U`

after each substep `m`

is

`Uᵐ⁺¹ = Uᵐ + Δt * (γᵐ * Gᵐ + ζᵐ * Gᵐ⁻¹)`

where `Uᵐ`

is the state at the $m$-th substep, `Gᵐ`

is the tendency at the $m$-th substep, `Gᵐ⁻¹`

is the tendency at the previous substep, and constants $γ¹ = 8/15$, $γ² = 5/12$, $γ³ = 3/4$, $ζ¹ = 0$, $ζ² = -17/60$, $ζ³ = -5/12$.

The state at the first substep is taken to be the one that corresponds to the $n$-th timestep, `U¹ = Uⁿ`

, and the state after the third substep is then the state at the `Uⁿ⁺¹ = U⁴`

.

`Oceananigans.TimeSteppers.time_step!`

— Method`time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false, compute_tendencies=true)`

Step forward `model`

one time step `Δt`

with a 2nd-order Adams-Bashforth method and pressure-correction substep. Setting `euler=true`

will take a forward Euler time step. Setting `compute_tendencies=false`

will not calculate new tendencies

`Oceananigans.TimeSteppers.time_step!`

— Method`time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt)`

Step forward `model`

one time step `Δt`

with a 3rd-order Runge-Kutta method. The 3rd-order Runge-Kutta method takes three intermediate substep stages to achieve a single timestep. A pressure correction step is applied at each intermediate stage.

## Turbulence closures

`Oceananigans.TurbulenceClosures.AnisotropicMinimumDissipation`

— Type`AnisotropicMinimumDissipation{FT} <: AbstractTurbulenceClosure`

Parameters for the "anisotropic minimum dissipation" turbulence closure for large eddy simulation proposed originally by Rozema *et al.* (2015) and Abkar *et al.* (2016), then modified by Verstappen (2018), and finally described and validated for by Vreugdenhil and Taylor (2018).

`Oceananigans.TurbulenceClosures.AnisotropicMinimumDissipation`

— Method```
AnisotropicMinimumDissipation([time_discretization = ExplicitTimeDiscretization, FT = Float64;]
C = 1/12, Cν = nothing, Cκ = nothing, Cb = nothing)
```

Return parameters of type `FT`

for the `AnisotropicMinimumDissipation`

turbulence closure.

**Arguments**

`time_discretization`

: Either`ExplicitTimeDiscretization()`

or`VerticallyImplicitTimeDiscretization()`

, which integrates the terms involving only $z$-derivatives in the viscous and diffusive fluxes with an implicit time discretization. Default`ExplicitTimeDiscretization()`

.`FT`

: Float type; default`Float64`

.

**Keyword arguments**

`C`

: Poincaré constant for both eddy viscosity and eddy diffusivities.`C`

is overridden for eddy viscosity or eddy diffusivity if`Cν`

or`Cκ`

are set, respecitvely.`Cν`

: Poincaré constant for momentum eddy viscosity.`Cκ`

: Poincaré constant for tracer eddy diffusivities. If one number or function, the same number or function is applied to all tracers. If a`NamedTuple`

, it must possess a field specifying the Poncaré constant for every tracer.`Cb`

: Buoyancy modification multiplier (`Cb = nothing`

turns it off,`Cb = 1`

was used by Abkar*et al.*(2016)).*Note*: that we*do not*subtract the horizontally-average component before computing this buoyancy modification term. This implementation differs from Abkar*et al.*(2016)'s proposal and the impact of this approximation has not been tested or validated.

By default: `C = Cν = Cκ = 1/12`

, which is appropriate for a finite-volume method employing a second-order advection scheme, and `Cb = nothing`

, which turns off the buoyancy modification term.

`Cν`

or `Cκ`

may be numbers, or functions of `x, y, z`

.

**Examples**

```
julia> using Oceananigans
julia> pretty_diffusive_closure = AnisotropicMinimumDissipation(C=1/2)
AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
Poincaré constant for momentum eddy viscosity Cν: 0.5
Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.5
Buoyancy modification multiplier Cb: nothing
```

```
julia> using Oceananigans
julia> const Δz = 0.5; # grid resolution at surface
julia> surface_enhanced_tracer_C(x, y, z) = 1/12 * (1 + exp((z + Δz/2) / 8Δz));
julia> fancy_closure = AnisotropicMinimumDissipation(Cκ=surface_enhanced_tracer_C)
AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: surface_enhanced_tracer_C
Buoyancy modification multiplier Cb: nothing
```

```
julia> using Oceananigans
julia> tracer_specific_closure = AnisotropicMinimumDissipation(Cκ=(c₁=1/12, c₂=1/6))
AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: (c₁ = 0.08333333333333333, c₂ = 0.16666666666666666)
Buoyancy modification multiplier Cb: nothing
```

**References**

Vreugdenhil C., and Taylor J. (2018), "Large-eddy simulations of stratified plane Couette flow using the anisotropic minimum-dissipation model", Physics of Fluids 30, 085104.

Verstappen, R. (2018), "How much eddy dissipation is needed to counterbalance the nonlinear production of small, unresolved scales in a large-eddy simulation of turbulence?", Computers & Fluids 176, pp. 276-284.

`Oceananigans.TurbulenceClosures.ConvectiveAdjustmentVerticalDiffusivity`

— Type```
ConvectiveAdjustmentVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(), FT=Float64;]
convective_κz = 0,
convective_νz = 0,
background_κz = 0,
background_νz = 0)
```

Return a convective adjustment vertical diffusivity closure that applies different values of diffusivity and/or viscosity depending whether the region is statically stable (positive or zero buoyancy gradient) or statically unstable (negative buoyancy gradient).

**Arguments**

`time_discretization`

: Either`ExplicitTimeDiscretization()`

or`VerticallyImplicitTimeDiscretization()`

; default`VerticallyImplicitTimeDiscretization()`

.`FT`

: Float type; default`Float64`

.

**Keyword arguments**

`convective_κz`

: Vertical tracer diffusivity in regions with negative (unstable) buoyancy gradients. Either a single number, function, array, field, or tuple of diffusivities for each tracer.`background_κz`

: Vertical tracer diffusivity in regions with zero or positive (stable) buoyancy gradients.`convective_νz`

: Vertical viscosity in regions with negative (unstable) buoyancy gradients. Either a number, function, array, or field.`background_κz`

: Vertical viscosity in regions with zero or positive (stable) buoyancy gradients.

**Example**

```
julia> using Oceananigans
julia> cavd = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1)
ConvectiveAdjustmentVerticalDiffusivity{VerticallyImplicitTimeDiscretization}(background_κz=0.0 convective_κz=1 background_νz=0.0 convective_νz=0.0)
```

`Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization`

— Type`struct ExplicitTimeDiscretization <: AbstractTimeDiscretization`

A fully-explicit time-discretization of a `TurbulenceClosure`

.

`Oceananigans.TurbulenceClosures.HorizontalDivergenceScalarDiffusivity`

— Type```
HorizontalDivergenceScalarDiffusivity([time_discretization=ExplicitTimeDiscretization(),
FT::DataType=Float64;]
kwargs...)
```

Shorthand for a `ScalarDiffusivity`

with `HorizontalDivergenceFormulation()`

. See `ScalarDiffusivity`

.

`Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity`

— Type```
HorizontalScalarDiffusivity([time_discretization=ExplicitTimeDiscretization(),
FT::DataType=Float64;]
kwargs...)
```

Shorthand for a `ScalarDiffusivity`

with `HorizontalFormulation()`

. See `ScalarDiffusivity`

.

`Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity`

— Method```
IsopycnalSkewSymmetricDiffusivity([time_disc=VerticallyImplicitTimeDiscretization(), FT=Float64;]
κ_skew = 0,
κ_symmetric = 0,
isopycnal_tensor = SmallSlopeIsopycnalTensor(),
slope_limiter = FluxTapering(1e-2))
```

Return parameters for an isopycnal skew-symmetric tracer diffusivity with skew diffusivity `κ_skew`

and symmetric diffusivity `κ_symmetric`

that uses an `isopycnal_tensor`

model for for calculating the isopycnal slopes, and (optionally) applying a `slope_limiter`

to the calculated isopycnal slope values.

Both `κ_skew`

and `κ_symmetric`

may be constants, arrays, fields, or functions of `(x, y, z, t)`

.

`Oceananigans.TurbulenceClosures.RiBasedVerticalDiffusivity`

— Type```
RiBasedVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(),
FT = Float64;]
Ri_dependent_tapering = HyperbolicTangentRiDependentTapering(),
ν₀ = 0.7,
κ₀ = 0.5,
κᶜᵃ = 1.7,
Cᵉⁿ = 0.1,
Cᵃᵛ = 0.6,
Ri₀ = 0.1,
Riᵟ = 0.4,
warning = true)
```

Return a closure that estimates the vertical viscosity and diffusivity from "convective adjustment" coefficients `ν₀`

and `κ₀`

multiplied by a decreasing function of the Richardson number, $Ri$.

**Arguments**

`time_discretization`

: Either`ExplicitTimeDiscretization()`

or`VerticallyImplicitTimeDiscretization()`

, which integrates the terms involving only $z$-derivatives in the viscous and diffusive fluxes with an implicit time discretization. Default`VerticallyImplicitTimeDiscretization()`

.`FT`

: Float type; default`Float64`

.

**Keyword arguments**

`Ri_dependent_tapering`

: The $Ri$-dependent tapering. Options are:`PiecewiseLinearRiDependentTapering()`

,`HyperbolicTangentRiDependentTapering()`

(default), and`ExponentialRiDependentTapering()`

.`ν₀`

: Non-convective viscosity.`κ₀`

: Non-convective diffusivity for tracers.`κᶜᵃ`

: Convective adjustment diffusivity for tracers.`Cᵉⁿ`

: Entrainment coefficient for tracers.`Cᵃᵛ`

: Time-averaging coefficient for viscosity and diffusivity.`Ri₀`

: $Ri$ threshold for decreasing viscosity and diffusivity.`Riᵟ`

: $Ri$-width over which viscosity and diffusivity decreases to 0.

`Oceananigans.TurbulenceClosures.ScalarBiharmonicDiffusivity`

— Type```
ScalarBiharmonicDiffusivity(formulation = ThreeDimensionalFormulation(), FT = Float64;
ν = 0,
κ = 0,
discrete_form = false,
loc = (nothing, nothing, nothing),
parameters = nothing)
```

Return a scalar biharmonic diffusivity turbulence closure with viscosity coefficient `ν`

and tracer diffusivities `κ`

for each tracer field in `tracers`

. If a single `κ`

is provided, it is applied to all tracers. Otherwise `κ`

must be a `NamedTuple`

with values for every tracer individually.

**Arguments**

`formulation`

:`HorizontalFormulation()`

for diffusivity applied in the horizontal direction(s)`VerticalFormulation()`

for diffusivity applied in the vertical direction,`ThreeDimensionalFormulation()`

(default) for diffusivity applied isotropically to all directions

`FT`

: the float datatype (default:`Float64`

)

**Keyword arguments**

`ν`

: Viscosity.`Number`

,`AbstractArray`

,`Field`

, or`Function`

.`κ`

: Diffusivity.`Number`

, three-dimensional`AbstractArray`

,`Field`

,`Function`

, or`NamedTuple`

of diffusivities with entries for each tracer.`discrete_form`

:`Boolean`

; default:`False`

.

When prescribing the viscosities or diffusivities as functions, depending on the value of keyword argument `discrete_form`

, the constructor expects:

`discrete_form = false`

(default): functions of the grid's native coordinates and time, e.g.,`(x, y, z, t)`

for a`RectilinearGrid`

or`(λ, φ, z, t)`

for a`LatitudeLongitudeGrid`

.`discrete_form = true`

:- with
`loc = (nothing, nothing, nothing)`

(default): functions of`(i, j, k, grid, ℓx, ℓy, ℓz)`

with`ℓx`

,`ℓy`

and`ℓz`

either`Face()`

or`Center()`

. - with
`loc = (ℓx, ℓy, ℓz)`

with`ℓx`

,`ℓy`

and`ℓz`

either`Face()`

or`Center()`

: functions of`(i, j, k, grid)`

.

- with
`parameters`

:`NamedTuple`

with parameters used by the functions that compute viscosity and/or diffusivity; default:`nothing`

.

For examples see `ScalarDiffusivity`

.

`Oceananigans.TurbulenceClosures.ScalarBiharmonicDiffusivity`

— Type`struct ScalarBiharmonicDiffusivity{F, N, K} <: AbstractScalarBiharmonicDiffusivity{F}`

Holds viscosity and diffusivities for models with prescribed isotropic diffusivities.

`Oceananigans.TurbulenceClosures.ScalarDiffusivity`

— Type```
ScalarDiffusivity(time_discretization = ExplicitTimeDiscretization(),
formulation = ThreeDimensionalFormulation(), FT = Float64;
ν = 0,
κ = 0,
discrete_form = false,
loc = (nothing, nothing, nothing),
parameters = nothing)
```

Return `ScalarDiffusivity`

turbulence closure with viscosity `ν`

and tracer diffusivities `κ`

for each tracer field in `tracers`

. If a single `κ`

is provided, it is applied to all tracers. Otherwise `κ`

must be a `NamedTuple`

with values for every tracer individually.

**Arguments**

```
* `time_discretization`: either `ExplicitTimeDiscretization()` (default)
or `VerticallyImplicitTimeDiscretization()`.
* `formulation`:
- `HorizontalFormulation()` for diffusivity applied in the horizontal direction(s)
- `VerticalFormulation()` for diffusivity applied in the vertical direction,
- `ThreeDimensionalFormulation()` (default) for diffusivity applied isotropically to all directions
* `FT`: the float datatype (default: `Float64`)
```

**Keyword arguments**

```
* `ν`: Viscosity. `Number`, three-dimensional `AbstractArray`, `Field`, or `Function`.
* `κ`: Diffusivity. `Number`, `AbstractArray`, `Field`, `Function`, or
`NamedTuple` of diffusivities with entries for each tracer.
* `discrete_form`: `Boolean`; default: `False`.
```

When prescribing the viscosities or diffusivities as functions, depending on the value of keyword argument `discrete_form`

, the constructor expects:

```
* `discrete_form = false` (default): functions of the grid's native coordinates
and time, e.g., `(x, y, z, t)` for a `RectilinearGrid` or
`(λ, φ, z, t)` for a `LatitudeLongitudeGrid`.
* `discrete_form = true`:
- with `loc = (nothing, nothing, nothing)` (default):
functions of `(i, j, k, grid, ℓx, ℓy, ℓz)` with `ℓx`, `ℓy`
and `ℓz` either `Face()` or `Center()`.
- with `loc = (ℓx, ℓy, ℓz)` with `ℓx`, `ℓy`
and `ℓz` either `Face()` or `Center()`: functions of `(i, j, k, grid)`.
* `parameters`: `NamedTuple` with parameters used by the functions
that compute viscosity and/or diffusivity; default: `nothing`.
```

**Examples**

```
julia> using Oceananigans
julia> ScalarDiffusivity(ν = 1000, κ=2000)
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1000.0, κ=2000.0)
```

```
julia> const depth_scale = 100;
julia> @inline ν(x, y, z) = 1000 * exp(z / depth_scale)
ν (generic function with 1 method)
julia> ScalarDiffusivity(ν = ν)
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=ν (generic function with 1 method), κ=0.0)
```

```
julia> using Oceananigans.Grids: znode
julia> @inline function κ(i, j, k, grid, ℓx, ℓy, ℓz)
z = znode(i, j, k, grid, ℓx, ℓy, ℓz)
return 2000 * exp(z / depth_scale)
end
κ (generic function with 1 method)
julia> ScalarDiffusivity(κ = κ, discrete_form = true)
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0, κ=Oceananigans.TurbulenceClosures.DiscreteDiffusionFunction{Nothing, Nothing, Nothing, Nothing, typeof(κ)})
```

```
julia> @inline function another_κ(i, j, k, grid, p)
z = znode(i, j, k, grid)
return 2000 * exp(z / p.depth_scale)
end
another_κ (generic function with 1 method)
julia> ScalarDiffusivity(κ = another_κ, discrete_form = true, loc = (Center, Center, Face), parameters = (; depth_scale = 120.0))
ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0, κ=Oceananigans.TurbulenceClosures.DiscreteDiffusionFunction{Center, Center, Face, NamedTuple{(:depth_scale,), Tuple{Float64}}, typeof(another_κ)})
```

`Oceananigans.TurbulenceClosures.SmagorinskyLilly`

— Method`SmagorinskyLilly([time_discretization::TD = ExplicitTimeDiscretization(), FT=Float64;] C=0.16, Cb=1.0, Pr=1.0)`

Return a `SmagorinskyLilly`

type associated with the turbulence closure proposed by Lilly (1962), Smagorinsky (1958), Smagorinsky (1963), and Lilly (1966), which has an eddy viscosity of the form

`νₑ = (C * Δᶠ)² * √(2Σ²) * √(1 - Cb * N² / Σ²)`

and an eddy diffusivity of the form

`κₑ = νₑ / Pr`

where `Δᶠ`

is the filter width, `Σ² = ΣᵢⱼΣᵢⱼ`

is the double dot product of the strain tensor `Σᵢⱼ`

, `Pr`

is the turbulent Prandtl number, `N²`

is the total buoyancy gradient, and `Cb`

is a constant the multiplies the Richardson number modification to the eddy viscosity.

**Arguments**

`time_discretization`

: Either`ExplicitTimeDiscretization()`

or`VerticallyImplicitTimeDiscretization()`

, which integrates the terms involving only $z$-derivatives in the viscous and diffusive fluxes with an implicit time discretization. Default`ExplicitTimeDiscretization()`

.`FT`

: Float type; default`Float64`

.

**Keyword arguments**

`C`

: Smagorinsky constant. Default value is 0.16 as obtained by Lilly (1966).`Cb`

: Buoyancy term multipler based on Lilly (1962) (`Cb = 0`

turns it off,`Cb ≠ 0`

turns it on. Typically, and according to the original work by Lilly (1962),`Cb = 1 / Pr`

.)`Pr`

: Turbulent Prandtl numbers for each tracer. Either a constant applied to every tracer, or a`NamedTuple`

with fields for each tracer individually.

**References**

Smagorinsky, J. "On the numerical integration of the primitive equations of motion for baroclinic flow in a closed region." Monthly Weather Review (1958)

Lilly, D. K. "On the numerical simulation of buoyant convection." Tellus (1962)

Smagorinsky, J. "General circulation experiments with the primitive equations: I. The basic experiment." Monthly Weather Review (1963)

Lilly, D. K. "The representation of small-scale turbulence in numerical simulation experiments." NCAR Manuscript No. 281, 0, (1966)

`Oceananigans.TurbulenceClosures.TwoDimensionalLeith`

— Type```
TwoDimensionalLeith(FT=Float64;
C=0.3, C_Redi=1, C_GM=1,
isopycnal_model=SmallSlopeIsopycnalTensor())
```

Return a `TwoDimensionalLeith`

type associated with the turbulence closure proposed by Leith (1968) and Fox-Kemper and Menemenlis (2008) which has an eddy viscosity of the form

`νₑ = (C * Δᶠ)³ * √(|∇ₕ ζ|² + |∇ₕ ∂w/∂z|²)`

and an eddy diffusivity of the form...

where `Δᶠ`

is the filter width, `ζ = ∂v/∂x - ∂u/∂y`

is the vertical vorticity, and `C`

is a model constant.

**Keyword arguments**

`C`

: Model constant`C_Redi`

: Coefficient for down-gradient tracer diffusivity for each tracer. Either a constant applied to every tracer, or a`NamedTuple`

with fields for each tracer individually.`C_GM`

: Coefficient for down-gradient tracer diffusivity for each tracer. Either a constant applied to every tracer, or a`NamedTuple`

with fields for each tracer individually.

**References**

Leith, C. E. (1968). "Diffusion Approximation for Two‐Dimensional Turbulence", The Physics of Fluids 11, 671. doi: 10.1063/1.1691968

Fox‐Kemper, B., & D. Menemenlis (2008), "Can large eddy simulation techniques improve mesoscale rich ocean models?", in Ocean Modeling in an Eddying Regime, Geophys. Monogr. Ser., 177, pp. 319–337. doi: 10.1029/177GM19

`Oceananigans.TurbulenceClosures.VerticalScalarDiffusivity`

— Type```
VerticalScalarDiffusivity([time_discretization=ExplicitTimeDiscretization(),
FT::DataType=Float64;]
kwargs...)
```

Shorthand for a `ScalarDiffusivity`

with `VerticalFormulation()`

. See `ScalarDiffusivity`

.

`Oceananigans.TurbulenceClosures.VerticallyImplicitTimeDiscretization`

— Type`struct VerticallyImplicitTimeDiscretization <: AbstractTimeDiscretization`

A vertically-implicit time-discretization of a `TurbulenceClosure`

.

This implies that a flux divergence such as $𝛁 ⋅ 𝐪$ at the $n$-th timestep is time-discretized as

`[∇ ⋅ q]ⁿ = [explicit_flux_divergence]ⁿ + [∂z (κ ∂z c)]ⁿ⁺¹`

`Oceananigans.TurbulenceClosures.diffusivity`

— Function`diffusivity(closure, tracer_index, diffusivity_fields)`

Returns the scalar diffusivity associated with `closure`

and `tracer_index`

.

`Oceananigans.TurbulenceClosures.viscosity`

— Function`viscosity(closure, diffusivities)`

Returns the scalar viscosity associated with `closure`

.

## Utilities

`Oceananigans.Utils.AndSchedule`

— Method`AndSchedule(schedules...)`

Return a schedule that actuates when all `child_schedule`

s actuate.

`Oceananigans.Utils.IterationInterval`

— Method`IterationInterval(interval; offset=0)`

Return a callable `IterationInterval`

that "actuates" (schedules output or callback execution) whenever the model iteration (modified by `offset`

) is a multiple of `interval`

.

For example,

`IterationInterval(100)`

actuates at iterations`[100, 200, 300, ...]`

.`IterationInterval(100, offset=-1)`

actuates at iterations`[99, 199, 299, ...]`

.

`Oceananigans.Utils.KernelParameters`

— Method`KernelParameters(size, offsets)`

Return parameters for kernel launching and execution that define (i) a tuple that defines the `size`

of the kernel being launched and (ii) a tuple of `offsets`

that offset loop indices. For example, `offsets = (0, 0, 0)`

with `size = (N, N, N)`

means all indices loop from `1:N`

. If `offsets = (1, 1, 1)`

, then all indices loop from `2:N+1`

. And so on.

**Example**

```
size = (8, 6, 4)
offsets = (0, 1, 2)
kp = KernelParameters(size, offsets)
# Launch a kernel with indices that range from i=1:8, j=2:7, k=3:6,
# where i, j, k are the first, second, and third index, respectively:
launch!(arch, grid, kp, kernel!; kernel_args...)
```

See the documentation for `launch!`

.

`Oceananigans.Utils.MultiRegionObject`

— Method`MultiRegionObject(regional_objects::Tuple; devices)`

Return a MultiRegionObject

`Oceananigans.Utils.OrSchedule`

— Method`OrSchedule(schedules...)`

Return a schedule that actuates when any of the `child_schedule`

s actuates.

`Oceananigans.Utils.SpecifiedTimes`

— Method`SpecifiedTimes(times)`

Return a callable `TimeInterval`

that "actuates" (schedules output or callback execution) whenever the model's clock equals the specified values in `times`

. For example,

`SpecifiedTimes([1, 15.3])`

actuates when`model.clock.time`

is`1`

and`15.3`

.

The specified `times`

need not be ordered as the `SpecifiedTimes`

constructor will check and order them in ascending order if needed.

`Oceananigans.Utils.TimeInterval`

— Type`struct TimeInterval <: AbstractSchedule`

Callable `TimeInterval`

schedule for periodic output or diagnostic evaluation according to `model.clock.time`

.

`Oceananigans.Utils.TimeInterval`

— Method`TimeInterval(interval)`

Return a callable `TimeInterval`

that schedules periodic output or diagnostic evaluation on a `interval`

of simulation time, as kept by `model.clock`

.

`Oceananigans.Utils.WallTimeInterval`

— Method`WallTimeInterval(interval; start_time = time_ns() * 1e-9)`

Return a callable `WallTimeInterval`

that schedules periodic output or diagnostic evaluation on a `interval`

of "wall time" while a simulation runs, in units of seconds.

The "wall time" is the actual real world time in seconds, as kept by an actual or hypothetical clock hanging on your wall.

The keyword argument `start_time`

can be used to specify a starting wall time other than the moment `WallTimeInterval`

is constructed.

`Oceananigans.Utils.launch!`

— Method`launch!(arch, grid, layout, kernel!, args...; kwargs...)`

Launches `kernel!`

, with arguments `args`

and keyword arguments `kwargs`

, over the `dims`

of `grid`

on the architecture `arch`

. kernels run on the default stream

`Oceananigans.Utils.pretty_filesize`

— Function`pretty_filesize(s, suffix="B")`

Convert a floating point value `s`

representing a file size to a more human-friendly formatted string with one decimal places with a `suffix`

defaulting to "B". Depending on the value of `s`

the string will be formatted to show `s`

using an SI prefix from bytes, kiB (1024 bytes), MiB (1024² bytes), and so on up to YiB (1024⁸ bytes).

`Oceananigans.Utils.prettytime`

— Function`prettytime(t, longform=true)`

Convert a floating point value `t`

representing an amount of time in SI units of seconds to a human-friendly string with three decimal places. Depending on the value of `t`

the string will be formatted to show `t`

in nanoseconds (ns), microseconds (μs), milliseconds (ms), seconds, minutes, hours, or days.

With `longform=false`

, we use s, m, hrs, and d in place of seconds, minutes, and hours.

`Oceananigans.Utils.with_tracers`

— Method`with_tracers(tracer_names, initial_tuple, tracer_default)`

Create a tuple corresponding to the solution variables `u`

, `v`

, `w`

, and `tracer_names`

. `initial_tuple`

is a `NamedTuple`

that at least has fields `u`

, `v`

, and `w`

, and may have some fields corresponding to the names in `tracer_names`

. `tracer_default`

is a function that produces a default tuple value for each tracer if not included in `initial_tuple`

.

`Oceananigans.Utils.work_layout`

— Method`work_layout(grid, dims; include_right_boundaries=false, location=nothing)`

Returns the `workgroup`

and `worksize`

for launching a kernel over `dims`

on `grid`

. The `workgroup`

is a tuple specifying the threads per block in each dimension. The `worksize`

specifies the range of the loop in each dimension.

Specifying `include_right_boundaries=true`

will ensure the work layout includes the right face end points along bounded dimensions. This requires the field `location`

to be specified.

For more information, see: https://github.com/CliMA/Oceananigans.jl/pull/308

`Oceananigans.Utils.@apply_regionally`

— Macro`@apply_regionally expr`

Distributes locally the function calls in `expr`

ession

It calls `apply_regionally!`

when the functions do not return anything.

In case the function in `expr`

returns something, `@apply_regionally`

calls `construct_regionally`

.

## Units

`Oceananigans.Units.GiB`

— Constant`GiB`

A `Float64`

constant equal to 1024`MiB`

. Useful for increasing the clarity of scripts, e.g. `max_filesize = 50GiB`

.

`Oceananigans.Units.KiB`

— Constant`KiB`

A `Float64`

constant equal to 1024.0. Useful for increasing the clarity of scripts, e.g. `max_filesize = 250KiB`

.

`Oceananigans.Units.MiB`

— Constant`MiB`

A `Float64`

constant equal to 1024`KiB`

. Useful for increasing the clarity of scripts, e.g. `max_filesize = 100MiB`

.

`Oceananigans.Units.TiB`

— Constant`TiB`

A `Float64`

constant equal to 1024`GiB`

. Useful for increasing the clarity of scripts, e.g. `max_filesize = 2TiB`

.

`Oceananigans.Units.day`

— Constant`day`

A `Float64`

constant equal to 24`hours`

. Useful for increasing the clarity of scripts, e.g. `stop_time = 1day`

.

`Oceananigans.Units.days`

— Constant`days`

A `Float64`

constant equal to 24`hours`

. Useful for increasing the clarity of scripts, e.g. `stop_time = 7days`

.

`Oceananigans.Units.hour`

— Constant`hour`

A `Float64`

constant equal to 60`minutes`

. Useful for increasing the clarity of scripts, e.g. `Δt = 1hour`

.

`Oceananigans.Units.hours`

— Constant`hours`

A `Float64`

constant equal to 60`minutes`

. Useful for increasing the clarity of scripts, e.g. `Δt = 3hours`

.

`Oceananigans.Units.kilometer`

— Constant`kilometer`

A `Float64`

constant equal to 1000`meters`

. Useful for increasing the clarity of scripts, e.g. `Lx = 1kilometer`

.

`Oceananigans.Units.kilometers`

— Constant`kilometers`

A `Float64`

constant equal to 1000`meters`

. Useful for increasing the clarity of scripts, e.g. `Lx = 5000kilometers`

.

`Oceananigans.Units.meter`

— Constant`meter`

A `Float64`

constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Lx = 1meter`

.

`Oceananigans.Units.meters`

— Constant`meters`

A `Float64`

constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Lx = 50meters`

.

`Oceananigans.Units.minute`

— Constant`minute`

A `Float64`

constant equal to 60`seconds`

. Useful for increasing the clarity of scripts, e.g. `Δt = 1minute`

.

`Oceananigans.Units.minutes`

— Constant`minutes`

A `Float64`

constant equal to 60`seconds`

. Useful for increasing the clarity of scripts, e.g. `Δt = 15minutes`

.

`Oceananigans.Units.second`

— Constant`second`

A `Float64`

constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Δt = 1second`

.

`Oceananigans.Units.seconds`

— Constant`seconds`

A `Float64`

constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Δt = 7seconds`

.

`Oceananigans.Units.Time`

— Type`Time(t)`

Return a time "selector" at the continuous time `t`

for linearly interpolating `FieldTimeSeries`

.

**Examples**

```
# Interpolate `field_time_series` to `t=0.1`, returning `interpolated::Field`
interpolated = field_time_series[Time(0.1)]
# Interpolate `field_time_series` at `i, j, k` and `t=0.1`
interpolated_ijk = field_time_series[i, j, k, Time(0.1)]
```