Library

AdvectionTerm(model; velocities, location)

Return a KernelFunctionOperation that computes the advection term, defined as

ADV = uᵢ∂ⱼ(uᵢuⱼ)

By default, the buoyancy production will be calculated using the resolved velocities and users cab use the keyword velocities to modify that behavior:

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));

julia> model = NonhydrostaticModel(grid=grid);

julia> using Oceanostics.TKEBudgetTerms: BuoyancyProductionTerm

julia> using Oceanostics.TKEBudgetTerms: AdvectionTerm

julia> ADV = AdvectionTerm(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢ∂ⱼuⱼuᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("NamedTuple", "Centered")

BuoyancyProductionTerm(model; velocities, tracers, location)

Return a KernelFunctionOperation that computes the buoyancy production term, defined as

BP = uᵢbᵢ

where bᵢ is the component of the buoyancy acceleration in the i-th direction (which is zero for x and y, except when gravity_unit_vector isn't aligned with the grid's z-direction) and all three components of i=1,2,3 are added up.

By default, the buoyancy production will be calculated using the resolved velocities and tracers:

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));

julia> model = NonhydrostaticModel(grid=grid, buoyancy=BuoyancyTracer(), tracers=:b);

julia> using Oceanostics.TKEBudgetTerms: BuoyancyProductionTerm

julia> wb = BuoyancyProductionTerm(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢbᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("NamedTuple", "BuoyancyForce", "NamedTuple")

If we want to calculate only the turbulent buoyancy production rate, we can do so by passing turbulent perturbations to the velocities and/or tracers options):

julia> w′ = Field(model.velocities.w - Field(Average(model.velocities.w)));

julia> b′ = Field(model.tracers.b - Field(Average(model.tracers.b)));

julia> w′b′ = BuoyancyProductionTerm(model, velocities=(u=model.velocities.u, v=model.velocities.v, w=w′), tracers=(b=b′,))
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢbᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("NamedTuple", "BuoyancyForce", "NamedTuple")

IsotropicKineticEnergyDissipationRate(
    model;
    U,
    V,
    W,
    location
)

Calculate the Viscous Dissipation Rate, defined as

ε = 2 ν SᵢⱼSᵢⱼ,

where Sᵢⱼ is the strain rate tensor, for a fluid with an isotropic turbulence closure (i.e., a turbulence closure where ν (eddy or not) is the same for all directions.

KineticEnergy(model, u, v, w; location, kwargs...)

Calculate the kinetic energy of model manually specifying u, v and w.

KineticEnergy(model; kwargs...)

Calculate the kinetic energy of model.

KineticEnergyDissipationRate(model; U, V, W, location)

Calculate the Kinetic Energy Dissipation Rate, defined as

ε = ν (∂uᵢ/∂xⱼ) (∂uᵢ/∂xⱼ)
ε = ∂ⱼuᵢ ⋅ Fᵢⱼ

where ∂ⱼuᵢ is the velocity gradient tensor and Fᵢⱼ is the stress tensor.

KineticEnergyForcingTerm(model; location)

Return a KernelFunctionOperation that computes the forcing term of the KE prognostic equation:

    FORC = uᵢFᵤᵢ

where uᵢ are the velocity components and Fᵤᵢ is the forcing term(s) in the uᵢ prognostic equation (i.e. the forcing for uᵢ).

KineticEnergyStressTerm(model; location)

Return a KernelFunctionOperation that computes the diffusive term of the KE prognostic equation:

    DIFF = uᵢ∂ⱼτᵢⱼ

where uᵢ are the velocity components and τᵢⱼ is the diffusive flux of i momentum in the j-th direction.

KineticEnergyTendency(model; location)

Return a KernelFunctionOperation that computes the tendency uᵢGᵢ of the KE, excluding the nonhydrostatic pressure contribution:

KET = ½∂ₜuᵢ² = uᵢGᵢ - uᵢ∂ᵢpₙₕₛ

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1, 1, 1));

julia> model = NonhydrostaticModel(; grid);

julia> using Oceanostics.TKEBudgetTerms: KineticEnergyTendency

julia> ke_tendency = KineticEnergyTendency(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢGᵢᶜᶜᶜ (generic function with 1 method)
└── arguments: ("Centered", "Nothing", "Nothing", "Nothing", "BoundaryCondition", "BoundaryCondition", "BoundaryCondition", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Nothing", "Clock", "NamedTuple")

PressureRedistributionTerm(
    model;
    velocities,
    pressure,
    location
)

Return a KernelFunctionOperation that computes the pressure redistribution term:

PR = uᵢ∂ᵢp

where p is the pressure. By default p is taken to be the total pressure (nonhydrostatic + hydrostatic):

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1,1,1));

julia> model = NonhydrostaticModel(grid=grid);

julia> using Oceanostics.TKEBudgetTerms: PressureRedistributionTerm

julia> ∇u⃗p = PressureRedistributionTerm(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢ∂ᵢpᶜᶜᶜ (generic function with 1 method)
└── arguments: ("NamedTuple", "Field")

We can also pass velocities and pressure keywords to perform more specific calculations. The example below illustrates calculation of the nonhydrostatic contribution to the pressure redistrubution term:

julia> ∇u⃗pNHS = PressureRedistributionTerm(model, pressure=model.pressures.pNHS)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: uᵢ∂ᵢpᶜᶜᶜ (generic function with 1 method)
└── arguments: ("NamedTuple", "Field")

TurbulentKineticEnergy(model, u, v, w; U, V, W, location)

Calculate the turbulent kinetic energy of model manually specifying u, v, w and optionally background velocities U, V and W.

TurbulentKineticEnergy(model; kwargs...)

Calculate the turbulent kinetic energy of model.

XShearProductionRate(model, u, v, w, U, V, W; location)

Calculate the shear production rate in the model's x direction, considering velocities u, v, w and background (or average) velocities U, V and W.

XShearProductionRate(model; U, V, W, kwargs...)

Calculate the shear production rate in the model's x direction. At least one of the mean velocities U, V and W must be specified otherwise the output will be zero.

YShearProductionRate(model, u, v, w, U, V, W; location)

Calculate the shear production rate in the model's y direction, considering velocities u, v, w and background (or average) velocities U, V and W.

YShearProductionRate(model; U, V, W, kwargs...)

Calculate the shear production rate in the model's y direction. At least one of the mean velocities U, V and W must be specified otherwise the output will be zero.

ZShearProductionRate(model, u, v, w, U, V, W; location)

Calculate the shear production rate in the model's z direction, considering velocities u, v, w and background (or average) velocities U, V and W.

ZShearProductionRate(model; U, V, W, kwargs...)

Calculate the shear production rate in the model's z direction. At least one of the mean velocities U, V and W must be specified otherwise the output will be zero.

ImmersedTracerDiffusion(
    model,
    c,
    c_immersed_bc,
    closure,
    diffusivity_fields,
    val_tracer_index,
    clock,
    model_fields;
    location
)

Calculates the diffusion term due to the bathymetry term as

DIFF = ∂ⱼ 𝓆ᶜⱼ,

where 𝓆ᶜⱼ is the bathymetry-led diffusion tensor for tracer c, using the Oceananigans' kernel immersed_∇_dot_qᶜ.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(; grid, tracers=:a);

julia> DIFF = ImmersedTracerDiffusion(model, :a)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: immersed_∇_dot_qᶜ (generic function with 2 methods)
└── arguments: ("Field", "BoundaryCondition", "Nothing", "Nothing", "Val", "Clock", "NamedTuple")

TotalTracerDiffusion(
    model,
    c,
    c_immersed_bc,
    closure,
    diffusivity_fields,
    val_tracer_index,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculates the total diffusion term as

DIFF = ∂ⱼ qᶜⱼ + ∂ⱼ 𝓆ᶜⱼ,

c. The calculation is done using the Oceananigans' kernels ∇_dot_qᶜ and immersed_∇_dot_qᶜ. where qᶜⱼ is the interior diffusion tensor and 𝓆ᶜⱼ is the bathymetry-led diffusion tensor for tracer

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(; grid, tracers=:a);

julia> DIFF = TotalTracerDiffusion(model, :a)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: total_∇_dot_qᶜ (generic function with 1 method)
└── arguments: ("Field", "BoundaryCondition", "Nothing", "Nothing", "Val", "Clock", "NamedTuple", "Nothing")

TracerAdvection(model, u, v, w, c, advection; location)

Calculates the advection of the tracer c as

ADV = ∂ⱼ (uⱼ c)

using Oceananigans' kernel div_Uc.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(; grid, tracers=:a);

julia> ADV = TracerAdvection(model, :a)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: div_Uc (generic function with 10 methods)
└── arguments: ("Centered", "NamedTuple", "Field")

TracerDiffusion(
    model,
    val_tracer_index,
    c,
    closure,
    diffusivity_fields,
    clock,
    model_fields,
    buoyancy;
    location
)

Calculates the diffusion term (excluding anything due to the bathymetry) as

DIFF = ∂ⱼ qᶜⱼ,

where qᶜⱼ is the diffusion tensor for tracer c, using the Oceananigans' kernel ∇_dot_qᶜ.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(; grid, tracers=:a);

julia> DIFF = TracerDiffusion(model, :a)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: ∇_dot_qᶜ (generic function with 10 methods)
└── arguments: ("Nothing", "Nothing", "Val", "Field", "Clock", "NamedTuple", "Nothing")

TracerForcing(model, forcing, clock, model_fields; location)

Calculate the forcing term Fᶜ on the equation for tracer c for model.

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(; grid, tracers=:a);

julia> FORC = TracerForcing(model, :a)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: zeroforcing (generic function with 1 method)
└── arguments: ("Clock", "NamedTuple")

TracerVarianceDiffusiveTerm(model, tracer_name; location)

Return a KernelFunctionOperation that computes the diffusive term of the tracer variance prognostic equation using Oceananigans' diffusive tracer flux divergence kernel:

    DIFF = 2 c ∂ⱼFⱼ

where c is the tracer, and Fⱼ is the tracer's diffusive flux in the j-th direction.

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, tracers=:b, closure=SmagorinskyLilly())

DIFF = TracerVarianceDiffusiveTerm(model, :b)

TracerVarianceDissipationRate(
    model,
    tracer_name;
    tracer,
    location
)

Return a KernelFunctionOperation that computes the isotropic variance dissipation rate for tracer_name in model.tracers. The isotropic variance dissipation rate is defined as

    χ = 2 ∂ⱼc ⋅ Fⱼ

where Fⱼ is the diffusive flux of c in the j-th direction and ∂ⱼ is the gradient operator. χ is implemented in its conservative formulation based on the equation above.

Note that often χ is written as χ = 2κ (∇c ⋅ ∇c), which is the special case for Fickian diffusion (κ is the tracer diffusivity).

Here tracer_name is needed even when passing tracer in order to get the appropriate tracer_index. When passing tracer, this function should be used as

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid=grid, tracers=:b, closure=SmagorinskyLilly())

b̄ = Field(Average(model.tracers.b, dims=(1,2)))
b′ = model.tracers.b - b̄

χb = TracerVarianceDissipationRate(model, :b, tracer=b′)

TracerVarianceTendency(model, tracer_name; location)

Return a KernelFunctionOperation that computes the tracer variance tendency:

TEND = 2 c ∂ₜc

where c is the tracer and ∂ₜc is the tracer tendency (computed using Oceananigans' tracer tendency kernel).

julia> using Oceananigans

julia> grid = RectilinearGrid(size = (1, 1, 4), extent = (1, 1, 1));

julia> model = NonhydrostaticModel(; grid, tracers=:b);

julia> using Oceanostics.TracerVarianceBudgetTerms: TracerVarianceTendency

julia> χ = TracerVarianceTendency(model, :b)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: c∂ₜcᶜᶜᶜ (generic function with 1 method)
└── arguments: ("Val", "Val", "Centered", "Nothing", "BoundaryCondition", "Nothing", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Clock", "typeof(Oceananigans.Forcings.zeroforcing)")

struct BuoyancyAnomalyCriterion{FT} <: Oceanostics.FlowDiagnostics.AbstractAnomalyCriterion

Defines the mixed layer to be the depth at which the buoyancy is more than threshold greater than the surface buoyancy (but the pertubaton is usually negative).

When this model is used, the arguments buoyancy_formulation and C should be supplied where C should be the named tuple (; b), with b the buoyancy tracer.

struct DensityAnomalyCriterion{FT} <: Oceanostics.FlowDiagnostics.AbstractAnomalyCriterion

Defines the mixed layer to be the depth at which the density is more than threshold greater than the surface density.

When this model is used, the arguments buoyancy_formulation and C should be supplied where buoyancy_formulation should be the buoyancy model, and C should be a named tuple of (; T, S), (; T) or (; S) (the latter two if the buoyancy model specifies a constant salinity or temperature).

BottomCellValue(diagnostic)

Returns the value of the given diagnostic at the bottom, which can be either the bottom of the domain (lowest vertical level) or an immersed bottom.

DirectionalErtelPotentialVorticity(
    model,
    direction;
    tracer_name,
    loc
)

Calculate the contribution from a given direction to the Ertel Potential Vorticity basde on a model and a direction. The Ertel Potential Vorticity is defined as

EPV = ωₜₒₜ ⋅ ∇b

where ωₜₒₜ is the total (relative + planetary) vorticity vector, b is the buoyancy and ∇ is the gradient operator.

ErtelPotentialVorticity(model; tracer_name, loc)

Calculate the Ertel Potential Vorticty for model, where the characteristics of the Coriolis rotation are taken from model.coriolis. The Ertel Potential Vorticity is defined as

EPV = ωₜₒₜ ⋅ ∇b

where ωₜₒₜ is the total (relative + planetary) vorticity vector, b is the buoyancy and ∇ is the gradient operator.

julia> using Oceananigans

julia> grid = RectilinearGrid(topology = (Flat, Flat, Bounded), size = 4, extent = 1);

julia> N² = 1e-6;

julia> b_bcs = FieldBoundaryConditions(top=GradientBoundaryCondition(N²));

julia> model = NonhydrostaticModel(; grid, coriolis=FPlane(1e-4), buoyancy=BuoyancyTracer(), tracers=:b, boundary_conditions=(; b=b_bcs));

julia> stratification(z) = N² * z;

julia> set!(model, b=stratification)

julia> using Oceanostics: ErtelPotentialVorticity

julia> EPV = ErtelPotentialVorticity(model)
KernelFunctionOperation at (Face, Face, Face)
├── grid: 1×1×4 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: ertel_potential_vorticity_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Field", "Int64", "Int64", "Float64")

julia> interior(compute!(Field(EPV)))
1×1×5 view(::Array{Float64, 3}, 1:1, 1:1, 4:8) with eltype Float64:
[:, :, 1] =
 0.0

[:, :, 2] =
 1.0000000000000002e-10

[:, :, 3] =
 9.999999999999998e-11

[:, :, 4] =
 1.0000000000000002e-10

[:, :, 5] =
 1.0e-10

Note that EPV values are correctly calculated both in the interior and the boundaries. In the interior and top boundary, EPV = f×N² = 10⁻¹⁰, while EPV = 0 at the bottom boundary since ∂b/∂z is zero there.

MixedLayerDepth(grid, args; criterion)

Returns the mixed layer depth defined as the depth at which criterion is true.

Defaults to DensityAnomalyCriterion where the depth is that at which the density is some threshold (defaults to 0.125kg/m³) higher than the surface density.

When DensityAnomalyCriterion is used, the arguments buoyancy_formulation and C should be supplied where buoyancy_formulation should be the buoyancy model, and C should be a named tuple of (; T, S), (; T) or (; S) (the latter two if the buoyancy model specifies a constant salinity or temperature).

QVelocityGradientTensorInvariant(model; loc)

Calculate the value of the Q velocity gradient tensor invariant. This is usually just called Q and it is generally used for identifying and visualizing vortices in fluid flow.

The definition and nomenclature comes from the equation for the eigenvalues λ of the velocity gradient tensor ∂ⱼuᵢ:

    λ³ + P λ² + Q λ + T = 0

from where Q is defined as

    Q = ½ (ΩᵢⱼΩᵢⱼ - SᵢⱼSᵢⱼ)

and where Sᵢⱼ= ½(∂ⱼuᵢ + ∂ᵢuⱼ) and Ωᵢⱼ= ½(∂ⱼuᵢ - ∂ᵢuⱼ). More info about it can be found in doi:10.1063/1.5124245.

RichardsonNumber(model; loc)

Calculate the Richardson Number as

    Ri = (∂b/∂z) / (|∂u⃗ₕ/∂z|²)

where z is the true vertical direction (ie anti-parallel to gravity).

RossbyNumber(
    model;
    loc,
    add_background,
    dWdy_bg,
    dVdz_bg,
    dUdz_bg,
    dWdx_bg,
    dUdy_bg,
    dVdx_bg
)

Calculate the Rossby number using the vorticity in the rotation axis direction according to model.coriolis. Rossby number is defined as

    Ro = ωᶻ / f

where ωᶻ is the vorticity in the Coriolis axis of rotation and f is the Coriolis rotation frequency.

StrainRateTensorModulus(model; loc)

Calculate the modulus (absolute value) of the strain rate tensor S, which is defined as the symmetric part of the velocity gradient tensor:

    Sᵢⱼ = ½(∂ⱼuᵢ + ∂ᵢuⱼ)

Its modulus is then defined (using Einstein summation notation) as

    || Sᵢⱼ || = √(Sᵢⱼ Sᵢⱼ)

ThermalWindPotentialVorticity(model; tracer_name, loc)

Calculate the Potential Vorticty assuming thermal wind balance for model, where the characteristics of the Coriolis rotation are taken from model.coriolis. The Potential Vorticity in this case is defined as

    TWPV = (f + ωᶻ) ∂b/∂z - f ((∂U/∂z)² + (∂V/∂z)²)

where f is the Coriolis frequency, ωᶻ is the relative vorticity in the z direction, b is the buoyancy, and ∂U/∂z and ∂V/∂z comprise the thermal wind shear.

VorticityTensorModulus(model; loc)

Calculate the modulus (absolute value) of the vorticity tensor Ω, which is defined as the antisymmetric part of the velocity gradient tensor:

    Ωᵢⱼ = ½(∂ⱼuᵢ - ∂ᵢuⱼ)

Its modulus is then defined (using Einstein summation notation) as

    || Ωᵢⱼ || = √(Ωᵢⱼ Ωᵢⱼ)

PotentialEnergy(model; location, geopotential_height)

Return a KernelFunctionOperation to compute the PotentialEnergy per unit volume,

\[Eₚ = \frac{gρ}{ρ₀}z = -bz\]

at each grid location in model. PotentialEnergy is defined for both BuoyancyTracer and SeawaterBuoyancy. See the relevant Oceananigans.jl documentation on buoyancy models for more information about available options.

The optional keyword argument geopotential_height is only used if ones wishes to calculate Eₚ with a potential density referenced to geopotential_height, rather than in-situ density, when using a BoussinesqEquationOfState.

Example

Usage with a BuoyancyTracer buoyacny model

julia> using Oceananigans

julia> using Oceanostics.PotentialEnergyEquationTerms: PotentialEnergy

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> model = NonhydrostaticModel(; grid, buoyancy=BuoyancyTracer(), tracers=(:b,))
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: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: minus_bz_ccc (generic function with 3 methods)
└── arguments: ("Field",)

The default behaviour of PotentialEnergy uses the in-situ density in the calculation when the equation of state is a BoussinesqEquationOfState:

julia> using Oceananigans, SeawaterPolynomials.TEOS10

julia> using Oceanostics.PotentialEnergyEquationTerms: PotentialEnergy

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: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: minus_bz_ccc (generic function with 3 methods)
└── arguments: ("KernelFunctionOperation", "NamedTuple")

To use a reference density set a constant value for the keyword argument geopotential_height and pass this the function. For example,

julia> using Oceananigans, SeawaterPolynomials.TEOS10;

julia> using Oceanostics.PotentialEnergyEquationTerms: PotentialEnergy;

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded));

julia> tracers = (:T, :S);

julia> eos = TEOS10EquationOfState();

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos);

julia> model = NonhydrostaticModel(; grid, buoyancy, tracers);

julia> geopotential_height = 0; # density variable will be σ₀

julia> PotentialEnergy(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: minus_bz_ccc (generic function with 3 methods)
└── arguments: ("KernelFunctionOperation", "NamedTuple")