API Documentation

Julia package for finite volume GPU and CPU large eddy simulations (LES) of atmospheric flows. The abstractions, design, and finite volume engine are based on Oceananigans.

AnelasticEquations

Module implementing anelastic dynamics for atmosphere models.

The anelastic approximation filters acoustic waves by assuming density and pressure are small perturbations from a dry, hydrostatic, adiabatic reference state. The key constraint is that mass flux divergence vanishes: ∇⋅(ρᵣ u) = 0.

AnelasticDynamics(
    reference_state
) -> AnelasticDynamics{_A, Nothing} where _A

Return AnelasticDynamics representing incompressible fluid dynamics expanded about reference_state.

AbstractMicrophysicalState{FT}

Abstract supertype for microphysical state structs.

Microphysical states encapsulate the local microphysical variables (e.g., cloud liquid, rain, droplet number) needed to compute tendencies. This abstraction enables the same tendency functions to work for both grid-based LES and Lagrangian parcel models.

Concrete subtypes should be immutable structs containing the relevant mixing ratios and number concentrations for a given microphysics scheme.

For example, a warm-phase one-moment scheme might define a state with cloud liquid and rain mixing ratios (qᶜˡ, qʳ).

See also microphysical_state, microphysical_tendency.

struct AllSkyOptics <: Breeze.AtmosphereModels.AbstractOptics

Type representing full-spectrum all-sky (cloudy) radiation using RRTMGP gas and cloud optics, can be used as optics argument in RadiativeTransferModel.

All-sky radiation includes scattering by cloud liquid and ice particles, requiring cloud water path, cloud fraction, and effective radius inputs from the microphysics scheme.

AtmosphereModel(
    grid;
    clock,
    thermodynamic_constants,
    formulation,
    dynamics,
    velocities,
    moisture_density,
    tracers,
    coriolis,
    boundary_conditions,
    forcing,
    advection,
    momentum_advection,
    scalar_advection,
    closure,
    microphysics,
    timestepper,
    timestepper_kwargs,
    radiation
) -> AtmosphereModel{Dyn, Frm, Arc, Tst, Grd, Clock{Float64, Float64, Int64, Int64}, Thm, Mom, Moi, Mfr, Nothing, Tmp, Sol, Vel, Trc, Adv, Nothing, Frc, Nothing, @NamedTuple{}, Nothing, Nothing, Nothing} where {Dyn, Frm, Arc, Tst, Grd, Thm, Mom, Moi, Mfr, Tmp, Sol, Vel, Trc, Adv, Frc}

Return an AtmosphereModel that uses the anelastic approximation following Pauluis (2008).

Arguments

• The default dynamics is AnelasticDynamics.

• The default formulation is :LiquidIcePotentialTemperature.

• The default advection scheme is Centered(order=2) for both momentum and scalars. If a single advection is provided, it is used for both momentum and scalars.

• Alternatively, specific momentum_advection and scalar_advection schemes may be provided. scalar_advection may be a NamedTuple with a different scheme for each respective scalar, identified by name.

Example

julia> using Breeze

julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 2, 3));

julia> model = AtmosphereModel(grid)
AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=288.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme:
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵗ: Centered(order=2)
├── forcing: @NamedTuple{ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵗ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

References

Pauluis, O. (2008). Thermodynamic consistency of the anelastic approximation for a moist atmosphere. Journal of the Atmospheric Sciences 65, 2719–2729.

struct AtmosphereModelBuoyancy{D, F, T}

Wrapper struct for computing buoyancy for AtmosphereModel in the context of a turbulence closure. Used to interface with Oceananigans turbulence closures that require buoyancy gradients.

struct BackgroundAtmosphere{FT}

Constant (spatially uniform) volume mixing ratios (VMR) for radiatively active gases. All values are dimensionless molar fractions.

Fields

• Major atmospheric constituents: N₂, O₂, CO₂, CH₄, N₂O, CO, NO₂, O₃
• Halocarbons: CFC₁₁, CFC₁₂, CFC₂₂, CCl₄, CF₄
• Hydrofluorocarbons: HFC₁₂₅, HFC₁₃₄ₐ, HFC₁₄₃ₐ, HFC₂₃, HFC₃₂

Defaults are approximate modern atmospheric values for major gases; halocarbons default to zero.

Note: H₂O is computed from the model's prognostic moisture field, not specified here.

struct ClearSkyOptics <: Breeze.AtmosphereModels.AbstractOptics

Type representing full-spectrum clear-sky radiation using RRTMGP gas optics, can be used as optics argument in RadiativeTransferModel.

struct ConstantRadiusParticles{FT}

• radius::Any: Effective radius [μm]

Represents cloud particles with a constant effective radius in microns (μm).

struct GrayOptics <: Breeze.AtmosphereModels.AbstractOptics

Type representing gray atmosphere radiation (O'Gorman & Schneider 2008), can be used as optics argument in RadiativeTransferModel.

struct HorizontalSlowMode{D}

Wrapper type indicating that vertical "fast" terms should be excluded from tendencies.

When computing momentum tendencies with a HorizontalSlowMode-wrapped dynamics, the horizontal pressure gradient is computed normally, but the vertical pressure gradient and buoyancy return zero. These vertical fast terms are handled by the acoustic substep loop through perturbation variables $-ψ ∂ρ''/∂z - g ρ''$.

Including the full vertical PG and buoyancy in the slow tendency introduces a hydrostatic truncation error $O(Δz^2)$ that drives spurious acoustic modes. The horizontal PG does not suffer from this issue and can safely be included.

NothingMicrophysicalState{FT}

A microphysical state with no prognostic variables.

Used for Nothing microphysics and SaturationAdjustment schemes where cloud condensate is diagnosed from the thermodynamic state rather than being prognostic.

RadiativeTransferModel(
    grid::Oceananigans.Grids.AbstractGrid,
    optics,
    args...;
    kw...
)

Construct a RadiativeTransferModel on grid using the specified optics.

Valid optics types are:

• GrayOptics() - Gray atmosphere radiation (O'Gorman & Schneider 2008)
• ClearSkyOptics() - Full-spectrum clear-sky radiation using RRTMGP gas optics
• AllSkyOptics() - Full-spectrum all-sky (cloudy) radiation using RRTMGP gas and cloud optics

The constants argument provides physical constants for the radiative transfer solver.

Example

julia> using Breeze, Oceananigans.Units, RRTMGP, NCDatasets

julia> grid = RectilinearGrid(; size=16, x=0, y=45, z=(0, 10kilometers),
                              topology=(Flat, Flat, Bounded));

julia> RadiativeTransferModel(grid, GrayOptics(), ThermodynamicConstants();
                              surface_temperature = 300,
                              surface_albedo = 0.1)
RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
└── diffuse_surface_albedo: ConstantField(0.1)

julia> RadiativeTransferModel(grid, ClearSkyOptics(), ThermodynamicConstants();
                              surface_temperature = 300,
                              surface_albedo = 0.1,
                              background_atmosphere = BackgroundAtmosphere(CO₂ = 400e-6))
RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
└── diffuse_surface_albedo: ConstantField(0.1)

struct SlowTendencyMode{D}

Wrapper type indicating that only "slow" tendencies should be computed.

When computing momentum tendencies with a SlowTendencyMode-wrapped dynamics, the "fast" terms (pressure gradient and buoyancy) return zero. This is used for split-explicit time-stepping where fast terms are handled separately in an acoustic substep loop.

See also SplitExplicitTimeDiscretization.

WarmRainState{FT} <: AbstractMicrophysicalState{FT}

A simple microphysical state for warm-rain schemes with cloud liquid and rain.

Fields

• qᶜˡ::Any: Specific cloud liquid water content [kg/kg]

• qʳ::Any: Specific rain water content [kg/kg]

buoyancy_forceᶜᶜᶜ(i, j, k, grid, dynamics, temperature,
                  specific_moisture, microphysics, microphysical_fields, constants)

Compute the buoyancy force density ρ b at cell center (i, j, k).

For anelastic dynamics, this returns -g (ρ - ρᵣ) where ρᵣ is the reference density. For compressible dynamics, this returns -g ρ directly.

This function is used in the vertical momentum equation to compute the gravitational forcing term.

compute_forcing!(forcing)

Compute any fields or quantities needed by a forcing before it is applied. This function is extended by the Forcings module for forcing types that require pre-computation (e.g., SubsidenceForcing which computes horizontal averages).

compute_pressure_correction!(model, Δt)

Compute the pressure correction for the given model. Default: no-op. For anelastic dynamics, solves the pressure Poisson equation.

dynamics_density(dynamics)

Return the density field appropriate to the dynamical formulation.

For anelastic dynamics, returns the reference density (time-independent background state). For compressible dynamics, returns the prognostic density field.

dynamics_pressure(dynamics)

Return the pressure field appropriate to the dynamical formulation.

For anelastic dynamics, returns the reference pressure (hydrostatic background state). For compressible dynamics, returns the prognostic pressure field.

grid_microphysical_tendency(i, j, k, grid, microphysics, name, ρ, fields, 𝒰, constants, velocities)

Compute the tendency for microphysical variable name at grid point (i, j, k).

This is the grid-indexed interface used by the tendency kernels. The default implementation builds the microphysical state ℳ via microphysical_state and dispatches to the state-based microphysical_tendency.

Schemes that need full grid access (e.g., for non-local operations) can override this method directly without using microphysical_state.

Arguments

• velocities: NamedTuple of velocity components (; u, v, w) [m/s].

initial_aerosol_number(microphysics) -> Any

Return the total initial aerosol number concentration [m⁻³] for a microphysics scheme.

This is used by initialize_model_microphysical_fields! and parcel model construction to set a physically meaningful default for the prognostic aerosol number density ρnᵃ. The value is derived from the aerosol size distribution stored in the microphysics scheme, so it stays consistent with the activation parameters.

Returns 0 by default; extensions override this for schemes with prognostic aerosol.

liquid_ice_potential_temperature(model)

Return the liquid-ice potential temperature field for the given model.

liquid_ice_potential_temperature_density(model)

Return the liquid-ice potential temperature density field for the given model.

make_pressure_correction!(model, Δt)

Apply the pressure correction to the momentum fields. Default: no-op. For anelastic dynamics, projects momentum to enforce the divergence constraint.

materialize_atmosphere_model_boundary_conditions(boundary_conditions, grid, formulation,
                                                dynamics, microphysics, surface_pressure, thermodynamic_constants,
                                                microphysical_fields, specific_moisture, temperature)

Regularize boundary conditions for an AtmosphereModel. This function is extended by the BoundaryConditions module to provide atmosphere-specific boundary condition handling.

If formulation is :LiquidIcePotentialTemperature and ρe boundary conditions are provided, they are automatically converted to ρθ boundary conditions by wrapping flux BCs in EnergyFluxBoundaryCondition, which divides by the local mixture heat capacity.

The dynamics argument provides access to the reference state for boundary conditions that require it, such as VirtualPotentialTemperature diagnostics.

The microphysics argument specifies the microphysics scheme used to compute moisture fractions for mixture heat capacity and virtual potential temperature calculations.

The microphysical_fields, specific_moisture, and temperature arguments are pre-created fields used to construct the VirtualPotentialTemperature diagnostic for stability-dependent boundary conditions.

materialize_atmosphere_model_forcing(forcing, field, name, model_field_names, context)

Materialize a forcing for an AtmosphereModel field. This function is extended by the Forcings module to handle atmosphere-specific forcing types like subsidence and geostrophic forcings.

The context argument provides additional information needed for materialization, such as grid, reference state, and thermodynamic constants.

mean_pressure(dynamics)

Return the mean (background/reference) pressure field in Pa.

microphysical_state(microphysics, ρ, μ, 𝒰, velocities)

Build an AbstractMicrophysicalState (ℳ) from density-weighted prognostic microphysical variables μ, density ρ, and thermodynamic state 𝒰.

This is the primary interface that microphysics schemes must implement. It converts density-weighted prognostics to the scheme-specific AbstractMicrophysicalState type.

For non-equilibrium schemes, cloud condensate comes from μ (prognostic fields). For saturation adjustment schemes, cloud condensate comes from 𝒰.moisture_mass_fractions, while precipitation (rain, snow) still comes from μ.

Arguments

• microphysics: The microphysics scheme
• ρ: Local density (scalar)
• μ: NamedTuple of density-weighted prognostic variables (e.g., (ρqᶜˡ=..., ρqʳ=...))
• 𝒰: Thermodynamic state
• velocities: NamedTuple of velocity components (; u, v, w) [m/s].

Returns

An AbstractMicrophysicalState subtype containing the local specific microphysical variables.

See also microphysical_tendency, AbstractMicrophysicalState.

microphysical_tendency(microphysics, name, ρ, ℳ, 𝒰, constants)

Compute the tendency for microphysical variable name from the microphysical state ℳ and thermodynamic state 𝒰.

This is the state-based tendency interface that operates on scalar states without grid indexing. It works identically for grid-based LES and parcel models.

Arguments

• microphysics: The microphysics scheme
• name: Variable name as Val(:name) (e.g., Val(:ρqᶜˡ))
• ρ: Local density (scalar)
• ℳ: Microphysical state (e.g., WarmPhaseOneMomentState)
• 𝒰: Thermodynamic state
• constants: Thermodynamic constants

Returns

The tendency value (scalar, units depend on variable).

See also microphysical_state, AbstractMicrophysicalState.

precipitation_rate(model, phase=:liquid)

Return a KernelFunctionOperation representing the precipitation rate for the given phase.

The precipitation rate is the rate at which moisture is removed from the atmosphere by precipitation processes. For zero-moment schemes, this is computed from the remove_precipitation function applied to cloud condensate.

Arguments:

• model: An AtmosphereModel with a microphysics scheme
• phase: Either :liquid (rain) or :ice (snow). Default is :liquid.

Returns a Field or KernelFunctionOperation that can be computed and visualized. Specific microphysics schemes must extend this function.

pressure_anomaly(dynamics)

Return the pressure anomaly (deviation from mean) in Pa.

set_to_mean!(reference_state, model)

Recompute the reference pressure and density profiles from horizontally-averaged temperature and moisture mass fractions of the current model state.

This is useful as a simulation callback to keep the reference state close to the evolving mean state, minimizing buoyancy perturbations ρ - ρᵣ.

specific_humidity(model) -> Any

Return the specific humidity (vapor mass fraction) field for the given model.

For Nothing microphysics (no condensate), the vapor mass fraction equals the total specific moisture. For microphysics schemes with prognostic vapor (e.g., where qᵛ is tracked explicitly), this function returns the appropriate vapor field.

static_energy(model)

Return the specific static energy field for the given model.

static_energy_density(model)

Return the static energy density field for the given model.

For LiquidIcePotentialTemperatureFormulation, returns a Field with boundary conditions that convert potential temperature fluxes to energy fluxes. This allows users to use BoundaryConditionOperation to extract energy flux values from the model.

For StaticEnergyFormulation, returns the prognostic energy density field directly.

surface_precipitation_flux(model) -> Any

Return a 2D Field representing the flux of precipitating moisture at the bottom boundary.

The surface precipitation flux is wʳ * ρqʳ at the bottom face (k=1), representing the rate at which rain mass leaves the domain through the bottom boundary.

Units: kg/m²/s (positive = downward flux out of domain)

Arguments:

• model: An AtmosphereModel with a microphysics scheme

Returns a 2D Field that can be computed and visualized. Specific microphysics schemes must extend this function.

thermodynamic_density(formulation)

Return the thermodynamic density field for the given formulation.

thermodynamic_density_name(formulation)

Return the name of the thermodynamic density field (e.g., :ρθ, :ρe, :ρE). Accepts a Symbol, Val(Symbol), or formulation struct.

total_pressure(dynamics)

Return the total pressure (mean + anomaly) in Pa.

Update auxiliary microphysical fields at grid point (i, j, k).

This is the single interface function for updating all auxiliary (non-prognostic) microphysical fields. Microphysics schemes should extend this function.

The function receives:

• μ: NamedTuple of microphysical fields (mutated)
• i, j, k: Grid indices (after μ since this is a mutating function)
• microphysics: The microphysics scheme
• ℳ: The microphysical state at this point
• ρ: Local density
• 𝒰: Thermodynamic state
• constants: Thermodynamic constants

Why i, j, k is needed

Grid indices cannot be eliminated because:

1. Fields must be written at specific grid points
2. Some schemes need grid-dependent logic (e.g., k == 1 for bottom boundary conditions in sedimentation schemes)

What to implement

Schemes should write all auxiliary fields in one function. This includes:

• Specific moisture fractions (qᶜˡ, qʳ, etc.) from the microphysical state
• Derived quantities (qˡ = qᶜˡ + qʳ, qⁱ = qᶜⁱ + qˢ)
• Vapor mass fraction qᵛ from the thermodynamic state
• Terminal velocities for sedimentation

See WarmRainState implementation below for an example.

update_microphysical_fields!(
    μ,
    i,
    j,
    k,
    grid,
    microphysics::Nothing,
    ρ,
    𝒰,
    constants
)

Update all microphysical fields at grid point (i, j, k).

This orchestrating function:

1. Builds the microphysical state ℳ via microphysical_state
2. Calls update_microphysical_auxiliaries! to write auxiliary fields

Schemes should implement update_microphysical_auxiliaries!, not this function.

DewpointTemperature(
    model;
    tolerance,
    maxiter
) -> KernelFunctionOperation{Center, Center, Center, _A, _B, K, Tuple{}} where {_A, _B, K<:(Breeze.AtmosphereModels.Diagnostics.DewpointTemperatureKernelFunction{_A, _B, _C, _D, _E, _F, Float64} where {_A, _B, _C, _D, _E, _F})}

Return a KernelFunctionOperation representing the dewpoint temperature $T⁺$.

The dewpoint temperature is the temperature at which the air would become saturated at its current vapor pressure. It is computed by solving the implicit equation:

\[pᵛ⁺(T⁺) = pᵛ\]

using secant iteration, where $pᵛ$ is the actual vapor pressure and $pᵛ⁺$ is the saturation vapor pressure.

For saturated air, the dewpoint temperature equals the actual temperature.

The keyword arguments tolerance (default 1e-4) and maxiter (default 10) control the secant iteration convergence.

Example

using Breeze

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))
model = AtmosphereModel(grid; microphysics=SaturationAdjustment())
set!(model, θ=300, qᵗ=0.01)

T⁺ = DewpointTemperature(model)

# output
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: DewpointTemperatureKernelFunction
└── arguments: ()

The result may be wrapped in a Field to store the computed values:

T⁺_field = Field(T⁺)

# output
1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11
    └── max=289.062, min=287.475, mean=288.27

EquivalentPotentialTemperature(model, flavor=:specific)

Return a KernelFunctionOperation representing equivalent potential temperature $θᵉ$.

Equivalent potential temperature is conserved during moist adiabatic processes (including condensation and evaporation) and is useful for identifying air masses and tracking convective processes. Following Emanuel1994 equation 4.5.11:

\[θᵉ = T \left(\frac{p₀}{p}\right)^{Rᵈ/cᵖᵐ} \exp\left(\frac{ℒˡ qᵛ}{cᵖᵐ T}\right) ℋ^γ\]

where $ℒˡ$ is the latent heat of vaporization, $qᵛ$ is the vapor mass fraction, $ℋ$ is the relative humidity, and $γ = -Rᵛ qᵛ / cᵖᵐ$.

Arguments

• model: An AtmosphereModel instance.
• flavor: Either :specific (default) to return $θᵉ$, or :density to return $ρ θᵉ$.

Examples

using Breeze

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))
model = AtmosphereModel(grid)
set!(model, θ=300, qᵗ=0.01)

θᵉ = EquivalentPotentialTemperature(model)
Field(θᵉ)

# output
1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11
    └── max=326.162, min=325.849, mean=326.005

References

• Emanuel, K. A. (1994). Atmospheric Convection. Oxford University Press.

LiquidIcePotentialTemperature(model, flavor=:specific)

Return a KernelFunctionOperation representing liquid-ice potential temperature $θˡⁱ$.

Liquid-ice potential temperature is a conserved quantity under moist adiabatic processes that accounts for the latent heat associated with liquid water and ice:

\[θˡⁱ = θ \left(1 - \frac{ℒˡᵣ qˡ + ℒⁱᵣ qⁱ}{cᵖᵐ T}\right)\]

where $θ$ is the mixture potential temperature, $ℒˡᵣ$ and $ℒⁱᵣ$ are the reference latent heats for liquid and ice, and $qˡ$, $qⁱ$ are the liquid and ice mass fractions.

Arguments

• model: An AtmosphereModel instance.
• flavor: Either :specific (default) to return $θˡⁱ$, or :density to return $ρ θˡⁱ$.

Examples

using Breeze

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))
model = AtmosphereModel(grid)
set!(model, θ=300, qᵗ=0.01)

θˡⁱ = LiquidIcePotentialTemperature(model)
Field(θˡⁱ)

# output
1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11
    └── max=300.0, min=300.0, mean=300.0

PotentialTemperature(model, flavor=:specific)

Return a KernelFunctionOperation representing the (mixture) potential temperature $θ$.

The potential temperature is defined as the temperature a parcel would have if adiabatically brought to a reference pressure $p₀$:

\[θ = \frac{T}{Π}\]

where $T$ is temperature and $Π = (p/p₀)^{Rᵐ/cᵖᵐ}$ is the mixture Exner function, computed using the moist air gas constant $Rᵐ$ and heat capacity $cᵖᵐ$.

Arguments

• model: An AtmosphereModel instance.
• flavor: Either :specific (default) to return $θ$, or :density to return $ρ θ$.

Examples

using Breeze

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))
model = AtmosphereModel(grid)
set!(model, θ=300, qᵗ=0.01)

θ = PotentialTemperature(model)
Field(θ)

# output
1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11
    └── max=300.0, min=300.0, mean=300.0

SaturationSpecificHumidity(
    model
) -> KernelFunctionOperation{Center, Center, Center, _A, _B, K, Tuple{}} where {_A, _B, K<:Breeze.AtmosphereModels.Diagnostics.SaturationSpecificHumidityKernelFunction}
SaturationSpecificHumidity(
    model,
    flavor_symbol
) -> KernelFunctionOperation{Center, Center, Center, _A, _B, K, Tuple{}} where {_A, _B, K<:Breeze.AtmosphereModels.Diagnostics.SaturationSpecificHumidityKernelFunction}

Return a KernelFunctionOperation representing the specified flavor of saturation specific humidity $qᵛ⁺$.

Flavor options

• :prognostic

  Return the saturation specific humidity corresponding to the model's prognostic state. This is the same as the equilibrium saturation specific humidity for saturated conditions and a model that uses saturation adjustment microphysics.

• :equilibrium

  Return the saturation specific humidity in potentially-saturated conditions, using the model.specific_moisture. This is equivalent to the :total_moisture flavor under saturated conditions with no condensate; or in other words, if model.specific_moisture happens to be equal to the saturation specific humidity.

• :total_moisture

  Return saturation specific humidity in the case that the total specific moisture is equal to the saturation specific humidity and there is no condensate. This is useful for manufacturing perfectly saturated initial conditions.

StabilityEquivalentPotentialTemperature(model, flavor=:specific)

Return a KernelFunctionOperation representing stability-equivalent potential temperature $θᵇ$.

Stability-equivalent potential temperature is a moist-conservative variable suitable for computing the moist Brunt-Väisälä frequency. It follows from the derivation in the paper by Durran and Klemp (1982), who show that the moist Brunt-Väisälä frequency $Nᵐ$ is correctly expressed in terms of the vertical gradient of a moist-conservative variable.

The formulation is based on equation (17) by Durran and Klemp (1982):

\[θᵇ = θᵉ \left( \frac{T}{Tᵣ} \right)^{cˡ qˡ / cᵖᵐ}\]

where $θᵉ$ is the equivalent potential temperature, $T$ is temperature, $Tᵣ$ is the energy reference temperature, $cˡ$ is the heat capacity of liquid water, $qᵗ$ is the total moisture specific humidity, and $cᵖᵐ$ is the moist air heat capacity.

This quantity is conserved along moist adiabats and is appropriate for use in stability calculations in saturated atmospheres.

Arguments

• model: An AtmosphereModel instance.
• flavor: Either :specific (default) to return $θᵇ$, or :density to return $ρ θᵇ$.

Examples

using Breeze

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))
model = AtmosphereModel(grid)
set!(model, θ=300, qᵗ=0.01)

θᵇ = StabilityEquivalentPotentialTemperature(model)
Field(θᵇ)

# output
1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11
    └── max=326.162, min=325.849, mean=326.005

References

• Durran, D. R. and Klemp, J. B. (1982). On the effects of moisture on the Brunt-Väisälä frequency. Journal of the Atmospheric Sciences 39, 2152–2158.

StaticEnergy(model, flavor=:specific)

Return a KernelFunctionOperation representing moist static energy $e$.

Moist static energy is a conserved quantity in adiabatic, frictionless flow that combines sensible heat, gravitational potential energy, and latent heat:

\[e = cᵖᵐ T + g z - ℒˡᵣ qˡ - ℒⁱᵣ qⁱ\]

where $cᵖᵐ$ is the moist air heat capacity, $T$ is temperature, $g$ is gravitational acceleration, $z$ is height, and $ℒˡᵣ qˡ + ℒⁱᵣ qⁱ$ is the latent heat content of condensate.

This is the prognostic thermodynamic variable used in StaticEnergyThermodynamics.

Arguments

• model: An AtmosphereModel instance.
• flavor: Either :specific (default) to return $e$, or :density to return $ρ e$.

Examples

using Breeze

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))
model = AtmosphereModel(grid)
set!(model, θ=300)

e = StaticEnergy(model)
Field(e)

# output
1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11
    └── max=3.03055e5, min=3.02663e5, mean=3.02859e5

VirtualPotentialTemperature(model, flavor=:specific)

Return a KernelFunctionOperation representing virtual potential temperature $θᵛ$.

Virtual potential temperature is the temperature that dry air would need to have in order to have the same density as moist air at the same pressure. To define virtual potential temperature, we first note the definition of virtual temperature:

\[Tᵛ = T \left( 1 + δᵛ qᵛ - qˡ - qⁱ \right)\]

where $δᵛ ≡ Rᵛ / Rᵈ - 1$. This follows from the ideal gas law for a mixture, $p = ρ Rᵐ T$, the mixture gas constant $Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ = Rᵈ \left( 1 + δᵛ qᵛ - qˡ - qⁱ \right)$, and the definition of virtual temperature, $p = ρ Rᵈ Tᵛ$, which leads to

\[Tᵛ = T rac{Rᵐ}{Rᵈ} = T \left( 1 + δᵛ qᵛ - qˡ - qⁱ \right)\]

The virtual potential temperature is defined analogously,

\[θᵛ = θ \left( 1 + δᵛ qᵛ - qˡ - qⁱ \right)\]

where $θ$ is potential temperature. Note that $Rᵛ / Rᵈ ≈ 1.608$ for water vapor and a dry air mixture typical to Earth's atmosphere, and that $δᵛ ≈ 0.608$.

using Breeze

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))
model = AtmosphereModel(grid)
set!(model, θ=300, qᵗ=0.01)

θᵛ = VirtualPotentialTemperature(model)
Field(θᵛ)

# output
1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×14 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11
    └── max=301.824, min=301.824, mean=301.824

BulkDragFunction(; direction=nothing, coefficient=1e-3, gustiness=0, surface_temperature=nothing)

Create a bulk drag function for computing surface momentum fluxes using bulk aerodynamic formulas. The drag function computes a quadratic drag:

\[Jᵘ = - Cᴰ |U| ρu\]

where Cᴰ is the drag coefficient, |U| = √(u² + v² + gustiness²) is the wind speed (with gustiness to prevent division by zero), and ρu is the momentum density.

Keyword Arguments

• direction: The direction of the momentum component (XDirection() or YDirection()). If nothing, the direction is inferred from the field location during boundary condition regularization.
• coefficient: The drag coefficient (default: 1e-3). Can be a constant or a PolynomialCoefficient for wind and stability-dependent transfer coefficients.
• gustiness: Minimum wind speed to prevent singularities when winds are calm (default: 0)
• surface_temperature: Surface temperature, required when using PolynomialCoefficient with stability correction. Can be a Field, Function, or Number. (default: nothing)

BulkSensibleHeatFluxFunction(; coefficient, gustiness=0, surface_temperature)

A bulk sensible heat flux function. The flux is computed as:

\[J = - ρ₀ Cᵀ |U| Δϕ\]

where Cᵀ is the transfer coefficient, |U| is the wind speed, and Δϕ is the difference between the near-surface atmospheric value and the surface value of the thermodynamic variable appropriate to the formulation:

• For LiquidIcePotentialTemperatureFormulation: Δϕ = θ - θ₀ (potential temperature flux)
• For StaticEnergyFormulation: Δϕ = e - cᵖᵈ T₀ (static energy flux)

The formulation is set automatically during model construction based on the thermodynamic formulation.

Keyword Arguments

• coefficient: The sensible heat transfer coefficient.
• gustiness: Minimum wind speed to prevent singularities (default: 0).
• surface_temperature: The surface temperature. Can be a Field, a Function, or a Number. Functions are converted to Fields during model construction.

BulkVaporFluxFunction(; coefficient, gustiness=0, surface_temperature)

Create a bulk vapor flux function for computing surface moisture fluxes. The flux is computed as:

\[Jᵛ = - ρ₀ Cᵛ |U| (qᵗ - qᵛ₀)\]

where Cᵛ is the transfer coefficient, |U| is the wind speed, qᵗ is the atmospheric specific humidity, and qᵛ₀ is the saturation specific humidity at the surface.

Keyword Arguments

• coefficient: The vapor transfer coefficient.
• gustiness: Minimum wind speed to prevent singularities (default: 0).
• surface_temperature: The surface temperature. Can be a Field, a Function, or a Number. Used to compute saturation specific humidity at the surface.

EnergyFluxBoundaryConditionFunction

A wrapper for boundary conditions that converts energy flux to potential temperature flux.

When using LiquidIcePotentialTemperatureFormulation, the prognostic thermodynamic variable is ρθ (potential temperature density). This wrapper allows users to specify energy fluxes (e.g., sensible heat flux in W/m²) which are converted to potential temperature fluxes by dividing by the local mixture heat capacity cᵖᵐ.

The relationship is:

\[Jᶿ = 𝒬 / cᵖᵐ\]

where 𝒬 is the energy flux and Jᶿ is the potential temperature flux.

The mixture heat capacity is computed using moisture fractions from the microphysics scheme, which correctly accounts for liquid and ice condensate when present.

FittedStabilityFunction(scalar_roughness_length;
    richardson_number_mapping = RichardsonNumberMapping(...),
    stability_function_parameters = StabilityFunctionParameters(...))

Stability correction based on Monin-Obukhov similarity theory using the Li et al. (2010) analytical mapping from bulk Richardson number to the stability parameter $ζ = z/L$.

Uses Hogström (1996) integrated stability functions for unstable conditions and Beljaars & Holtslag (1991) for stable conditions.

Applies structurally correct (and different) corrections for momentum vs scalar transfer:

• Momentum: $Cᴰ = Cᴰ_N [α / (α - Ψᴰ)]²$
• Scalar: $Cᵀ = Cᵀ_N [α / (α - Ψᴰ)] [β_h / (β_h - Ψᵀ)]$

where $α = \ln(z/ℓ)$, $β_h = \ln(z/ℓ_h)$.

FittedStabilityFunction is callable: sf(Riᴮ, α, β) returns the momentum stability correction factor, and sf(Riᴮ, α, β, Val(:scalar)) returns the scalar correction factor.

Arguments

• scalar_roughness_length: Roughness length for heat/moisture $ℓ_h$ (m).

Keyword Arguments

• richardson_number_mapping: RichardsonNumberMapping coefficients (default: Li et al. 2010).
• stability_function_parameters: StabilityFunctionParameters (default: Hogström 1996 / Beljaars & Holtslag 1991).

References

• Li, Y., Gao, Z., Lenschow, D. H., & Chen, F. (2010). An improved approach for parameterizing surface-layer turbulent transfer coefficients in numerical models. Boundary-Layer Meteorology, 137, 153-165.
• Hogström, U. (1996). Review of some basic characteristics of the atmospheric surface layer. Boundary-Layer Meteorology, 78, 215-246.
• Beljaars, A. C. M., & Holtslag, A. A. M. (1991). Flux parameterization over land surfaces for atmospheric models. Journal of Applied Meteorology, 30, 327-341.

PolynomialCoefficient(;
    polynomial = nothing,
    roughness_length = 1.5e-4,
    stability_function = FittedStabilityFunction(roughness_length / 7.3),
    surface = PlanarLiquidSurface()
)

A bulk transfer coefficient that depends on wind speed and atmospheric stability, following Large and Yeager (2009).

The neutral transfer coefficient at 10 m follows the Large & Yeager (2009) form:

\[C^N_{10}(U_h) = (a_0 + a_1 U_h + a_2 / U_h) × 10^{-3}\]

where $U_h$ is the wind speed at measurement height $h$.

The coefficient is adjusted for measurement height using logarithmic profile theory, and stability correction is applied based on the bulk Richardson number.

When polynomial is nothing, the appropriate Large & Yeager (2009) polynomial will be automatically selected based on the boundary condition type:

• BulkDrag: default_neutral_drag_polynomial = (0.142, 0.076, 2.7) for momentum
• BulkSensibleHeatFlux: default_neutral_sensible_heat_polynomial = (0.128, 0.068, 2.43) for sensible heat
• BulkVaporFlux: default_neutral_latent_heat_polynomial = (0.120, 0.070, 2.55) for latent heat

Keyword Arguments

• polynomial: Tuple (a₀, a₁, a₂) for the polynomial. If nothing, the polynomial is automatically selected by the boundary condition constructor.
• roughness_length: Surface roughness ℓ in meters (default: 1.5e-4, typical for ocean)
• minimum_wind_speed: Minimum wind speed to avoid singularity in a₂/U term (default: 0.1 m/s)
• stability_function: Stability correction strategy. Default is FittedStabilityFunction using Li et al. (2010) $Riᴮ → ζ$ mapping with Hogström (1996) / Beljaars & Holtslag (1991) MOST stability functions. The scalar roughness length defaults to roughness_length / 7.3 (typical ocean value). Use nothing to disable stability correction.
• surface: Surface type for computing saturation specific humidity in the stability correction. Default is PlanarLiquidSurface(). Use PlanarIceSurface() for ice surfaces.

The measurement height is automatically determined from the grid as the height of the first cell center above the surface.

Examples

using Breeze.BoundaryConditions: PolynomialCoefficient

# Polynomial coefficient with default settings
coef = PolynomialCoefficient()

# output
PolynomialCoefficient{Float64}
├── polynomial: nothing
├── roughness_length: 0.00015 m
├── minimum_wind_speed: 0.1 m/s
├── surface: PlanarLiquidSurface
└── stability_function: FittedStabilityFunction (Li et al. 2010)

using Breeze.BoundaryConditions: PolynomialCoefficient

# With explicit polynomial
coef = PolynomialCoefficient(polynomial = (0.142, 0.076, 2.7))

# output
PolynomialCoefficient{Float64}
├── polynomial: (0.142, 0.076, 2.7)
├── roughness_length: 0.00015 m
├── minimum_wind_speed: 0.1 m/s
├── surface: PlanarLiquidSurface
└── stability_function: FittedStabilityFunction (Li et al. 2010)

using Breeze.BoundaryConditions: PolynomialCoefficient

# No stability correction
coef = PolynomialCoefficient(stability_function = nothing)

# output
PolynomialCoefficient{Float64}
├── polynomial: nothing
├── roughness_length: 0.00015 m
├── minimum_wind_speed: 0.1 m/s
├── surface: PlanarLiquidSurface
└── stability_function: Nothing

References

• Large, W., & Yeager, S. G. (2009). The global climatology of an interannually varying air–sea flux data set. Climate dynamics, 33(2), 341-364.
• Li, Y., Gao, Z., Lenschow, D. H., & Chen, F. (2010). An improved approach for parameterizing surface-layer turbulent transfer coefficients in numerical models. Boundary-Layer Meteorology, 137, 153-165.

Evaluate the bulk transfer coefficient for given conditions.

For a materialized PolynomialCoefficient (with virtual_potential_temperature, surface_pressure, and thermodynamic_constants filled in during model construction), the stability correction is computed internally from the stored fields.

Arguments

• i, j: Grid indices
• grid: The grid
• U: Wind speed (m/s)
• T₀: Surface temperature (K) at location (i, j)

Returns the transfer coefficient (dimensionless).

RichardsonNumberMapping(FT = Oceananigans.defaults.FloatType; kwargs...)

Regression coefficients for the non-iterative mapping from bulk Richardson number $Riᴮ$ to the Monin-Obukhov stability parameter $ζ = z/L$, following Li et al. (2010).

The superscripts u, w, s denote unstable, weakly stable, and strongly stable regimes respectively. Subscript indices follow the original paper.

Three regimes:

• Unstable ($Riᴮ <$ stable_unstable_transition): Eq. 12
• Weakly stable (stable_unstable_transition $≤ Riᴮ ≤$ strongly_stable_transition): Eq. 14
• Strongly stable ($Riᴮ >$ strongly_stable_transition): Eq. 16

References

• Li, Y., Gao, Z., Lenschow, D. H., & Chen, F. (2010). An improved approach for parameterizing surface-layer turbulent transfer coefficients in numerical models. Boundary-Layer Meteorology, 137, 153-165.

StabilityFunctionParameters(FT = Oceananigans.defaults.FloatType;
    γᴰ = 19.3, γᵀ = 11.6, a = 1, b = 2/3, c = 5, d = 0.35)

Parameters for the integrated Monin-Obukhov stability functions $Ψ^D(ζ)$ and $Ψ^T(ζ)$.

Note: we use superscript D (drag/momentum) and T (temperature/scalar) to match the transfer coefficient notation $Cᴰ$, $Cᵀ$ established in notation.md. In the literature these are commonly written $Ψ_m$ and $Ψ_h$.

For unstable conditions ($ζ < 0$), uses Hogström (1996):

• $φ^D = (1 - γ^D ζ)^{-1/4}$
• $φ^T = 0.95(1 - γ^T ζ)^{-1/2}$

For stable conditions ($ζ ≥ 0$), uses Beljaars & Holtslag (1991):

• $Ψ^D = -(a ζ + b (ζ - c/d) e^{-dζ} + bc/d)$
• $Ψ^T = -((1 + 2aζ/3)^{3/2} + b (ζ - c/d) e^{-dζ} + bc/d - 1)$

References

• Hogström, U. (1996). Review of some basic characteristics of the atmospheric surface layer. Boundary-Layer Meteorology, 78, 215-246.
• Beljaars, A. C. M., & Holtslag, A. A. M. (1991). Flux parameterization over land surfaces for atmospheric models. Journal of Applied Meteorology, 30, 327-341.

ThetaFluxBoundaryConditionFunction

A wrapper for boundary conditions that converts potential temperature flux to energy flux.

When building a diagnostic energy_density field from a PotentialTemperatureFormulation, the boundary conditions on ρθ (potential temperature density) must be converted to energy flux boundary conditions by multiplying by the local mixture heat capacity cᵖᵐ.

The relationship is:

\[𝒬 = Jᶿ × cᵖᵐ\]

where 𝒬 is the energy flux and Jᶿ is the potential temperature flux.

BulkDrag(; direction=nothing, coefficient=1e-3, gustiness=0, surface_temperature=nothing)

Create a FluxBoundaryCondition for surface momentum drag.

See BulkDragFunction for details.

Examples

using Breeze

drag = BulkDrag(coefficient=1e-3, gustiness=0.1)

# output
FluxBoundaryCondition: BulkDragFunction(direction=Nothing, coefficient=0.001, gustiness=0.1)

Or with explicit direction, e.g., XDirection() for u:

using Oceananigans.Grids: XDirection

u_drag = BulkDrag(direction=XDirection(), coefficient=1e-3)
ρu_bcs = FieldBoundaryConditions(bottom=u_drag)

# output
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: FluxBoundaryCondition: BulkDragFunction(direction=XDirection(), coefficient=0.001, gustiness=0)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

and similarly for YDirection for v.

BulkSensibleHeatFlux(; coefficient, gustiness=0, surface_temperature)

Create a FluxBoundaryCondition for surface sensible heat flux.

The bulk formula computes J = -ρ₀ Cᵀ |U| Δϕ, where Δϕ depends on the thermodynamic formulation: Δθ for potential temperature or Δe for static energy. The formulation is set automatically during model construction.

See BulkSensibleHeatFluxFunction for details.

Example

using Breeze

T₀(x, y) = 290 + 2 * sign(cos(2π * x / 20e3))

ρe_bc = BulkSensibleHeatFlux(coefficient = 1e-3,
                             gustiness = 0.1,
                             surface_temperature = T₀)

# output
FluxBoundaryCondition: BulkSensibleHeatFluxFunction(coefficient=0.001, gustiness=0.1)

BulkVaporFlux(; coefficient, surface_temperature, gustiness=0)

Create a FluxBoundaryCondition for surface moisture flux.

The saturation specific humidity at the surface is automatically computed from surface_temperature.

See BulkVaporFluxFunction for details.

Example

using Breeze

T₀(x, y) = 290 + 2 * sign(cos(2π * x / 20e3))

ρqᵗ_bc = BulkVaporFlux(coefficient = 1e-3,
                       gustiness = 0.1,
                       surface_temperature = T₀)

# output
FluxBoundaryCondition: BulkVaporFluxFunction(coefficient=0.001, gustiness=0.1)

EnergyFluxBoundaryCondition(flux)

Create a boundary condition that wraps an energy flux and converts it to a potential temperature flux for use with LiquidIcePotentialTemperatureFormulation.

The energy flux is divided by the local mixture heat capacity cᵖᵐ to obtain the potential temperature flux: Jᶿ = 𝒬 / cᵖᵐ.

ThetaFluxBoundaryCondition(flux)

Create a boundary condition that wraps a potential temperature flux and converts it to an energy flux for use with diagnostic energy density fields.

The potential temperature flux is multiplied by the local mixture heat capacity cᵖᵐ to obtain the energy flux: 𝒬 = Jᶿ × cᵖᵐ.

cos_solar_zenith_angle(
    i,
    j,
    grid::RectilinearGrid{<:Any, <:Flat, <:Flat, <:Bounded},
    datetime::Dates.DateTime
) -> Any

Compute the cosine of the solar zenith angle for the grid's location.

For single-column grids with Flat horizontal topology, extracts latitude from the y-coordinate and longitude from the x-coordinate.

cos_solar_zenith_angle(
    datetime::Dates.DateTime,
    longitude,
    latitude
) -> Any

Compute the cosine of the solar zenith angle for a given datetime and location.

The solar zenith angle $θ_z$ satisfies:

\[\cos(θ_z) = \sin(φ) \sin(δ) + \cos(φ) \cos(δ) \cos(ω)\]

where:

• $φ$ is the latitude
• $δ$ is the solar declination
• $ω$ is the hour angle

Arguments

• datetime: UTC datetime
• latitude: latitude in degrees (positive North)
• longitude: longitude in degrees (positive East)

Returns

A value between -1 and 1. Negative values indicate the sun is below the horizon.

day_of_year(dt::Dates.DateTime) -> Int64

Return the day of year (1-365/366) for a given DateTime.

equation_of_time(day_of_year) -> Any

Compute the equation of time (in minutes) for a given day of year.

This accounts for the difference between mean solar time and apparent solar time due to the eccentricity of Earth's orbit and the obliquity of the ecliptic.

Uses the approximation by Spencer (1971); see solar_declination.

References

• Spencer, J. W. (1971) Fourier series representation of the position of the sun. Search, 2, 162-172.

hour_angle(datetime::Dates.DateTime, longitude) -> Any

Compute the hour angle (in radians) for a given datetime and longitude.

The hour angle $ω$ is zero at solar noon and increases by 15° per hour (Earth rotates 360°/24h = 15°/h).

Arguments

• datetime: UTC datetime
• longitude: longitude in degrees (positive East)

solar_declination(day_of_year) -> Any

Compute the solar declination angle (in radians) for a given day of year.

Uses the approximation by Spencer (1971):

\[δ = 0.006918 - 0.399912 \cos(γ) + 0.070257 \sin(γ) - 0.006758 \cos(2γ) + 0.000907 \sin(2γ) - 0.002697 \cos(3γ) + 0.00148 \sin(3γ)\]

where $γ = 2π (d - 1) / 365$ is the fractional year in radians and $d$ is the day of year.

References

• Spencer, J. W. (1971) Fourier series representation of the position of the sun. Search, 2, 162-172.

CompressibleEquations

Module implementing fully compressible dynamics for atmosphere models.

The compressible formulation directly time-steps density as a prognostic variable and computes pressure from the ideal gas law. This formulation does not filter acoustic waves, so explicit time-stepping with small time steps (or acoustic substepping) is required.

The fully compressible Euler equations in conservation form are:

\[\begin{aligned} &\text{Mass:} && \partial_t \rho + \boldsymbol{\nabla \cdot} (\rho \boldsymbol{u}) = 0 \\ &\text{Momentum:} && \partial_t (\rho \boldsymbol{u}) + \boldsymbol{\nabla \cdot} (\rho \boldsymbol{u} \boldsymbol{u}) + \boldsymbol{\nabla} p = -\rho g \hat{\boldsymbol{z}} + \rho \boldsymbol{f} + \boldsymbol{\nabla \cdot \mathcal{T}} \end{aligned}\]

Pressure is computed from the ideal gas law:

\[p = \rho R^m T\]

where $R^m$ is the mixture gas constant.

AcousticSubstepper

Storage and parameters for acoustic substepping using the Exner pressure formulation, following CM1's sound.F.

The acoustic loop uses velocity (u, v, w) and Exner pressure perturbation (π') as prognostic variables, forming a stable 2-variable system that avoids the density-buoyancy instability of the (ρu, ρ, ρθ) formulation.

The forward-backward scheme updates:

1. Forward: Velocity from Exner pressure gradient: $u += Δτ (u_{ten} - cᵖ θᵥ ∂π'_d/∂x)$
2. Backward: Exner pressure from velocity divergence: $π' += Δτ (π_{ten} - S ∇·u) + w_{terms}$
3. Implicit: Vertically implicit w-π' coupling (tridiagonal solve)
4. Filtering: $π̃' = π' + ϰᵈⁱ (π' - π'_{old})$

Fields

• substeps: Number of acoustic substeps for the full time step
• forward_weight: Off-centering parameter ω (0.6 = CM1 default)
• divergence_damping_coefficient: Forward-extrapolation filter ϰᵈⁱ for π' (default 0.10)
• acoustic_damping_coefficient: Klemp 2018 ϰᵃᶜ for velocity damping
• virtual_potential_temperature: Stage-frozen θᵥ (CenterField)
• acoustic_compression: Coefficient (γ-1)π₀ converting ∇·u to ∂π'/∂t (CenterField)
• reference_exner_function: Reference π₀ = (p_ref/pˢᵗ)^(R/cᵖ) (CenterField)
• exner_perturbation: Current Exner pressure perturbation π' = π - π₀ (CenterField)
• previous_exner_perturbation: Previous-substep π' for divergence damping (CenterField)
• filtered_exner_perturbation: Filtered π̃' used in PGF (CenterField)
• stage_thermodynamic_density: Stage-frozen ρθ (CenterField)
• averaged_velocities: Time-averaged velocities for scalar advection
• slow_tendencies: Frozen slow tendencies (velocity, exner_pressure). Momentum tendencies are stored in the outer timestepper's Gⁿ fields; density and thermodynamic density tendencies are also read directly from Gⁿ.
• vertical_solver: BatchedTridiagonalSolver for implicit w-π' coupling
• rhs: Right-hand side storage for tridiagonal solve

AcousticSubstepper(
    grid,
    split_explicit::SplitExplicitTimeDiscretization
) -> AcousticSubstepper{_A, _B, CF, AV, ST, TS} where {_A, _B, CF<:(Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), AV<:(NamedTuple{(:u, :v, :w), <:Tuple{Field{Face, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Face, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Face, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}}), ST<:(NamedTuple{(:velocity, :exner_pressure), <:Tuple{NamedTuple{(:u, :v, :w), <:Tuple{Field{Face, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Face, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Face, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}}, Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}}), TS<:(Oceananigans.Solvers.BatchedTridiagonalSolver{_A, _B, _C, _D, G, Nothing, Oceananigans.Grids.ZDirection} where {_A, _B, _C, _D, G<:Oceananigans.Grids.AbstractGrid})}

Construct an AcousticSubstepper using the Exner pressure formulation.

struct CompressibleDynamics{TD, D, P, FT, RS}

Fully compressible dynamics with prognostic density and diagnostic pressure.

Fields

• density: Prognostic density field ρ
• pressure: Diagnostic pressure field p = ρ Rᵐ T
• standard_pressure: Reference pressure pˢᵗ for potential temperature (default 10⁵ Pa)
• surface_pressure: Mean pressure at the bottom of the atmosphere p₀
• time_discretization: Time discretization scheme (SplitExplicitTimeDiscretization or ExplicitTimeStepping)
• reference_state: Fixed hydrostatically-balanced reference state for base-state pressure correction (nothing or ExnerReferenceState)

The time_discretization determines how tendencies are computed and which time-stepper is used:

• SplitExplicitTimeDiscretization: Acoustic substepping with separate slow/fast tendencies
• ExplicitTimeStepping: All tendencies computed together (small Δt required)

CompressibleDynamics(
;
    ...
) -> CompressibleDynamics{ExplicitTimeStepping, Nothing, Nothing, Float64, Nothing}
CompressibleDynamics(
    time_discretization;
    standard_pressure,
    surface_pressure,
    reference_potential_temperature
) -> CompressibleDynamics{_A, Nothing, Nothing, Float64, Nothing} where _A

Construct CompressibleDynamics. The density and pressure fields are materialized later in the model constructor.

Positional Arguments

• time_discretization: Time discretization scheme. Default: ExplicitTimeStepping. Use SplitExplicitTimeDiscretization for acoustic substepping.

Keyword Arguments

• standard_pressure: Reference pressure for potential temperature (default: 10⁵ Pa)
• surface_pressure: Mean surface pressure (default: 101325.0 Pa)
• reference_potential_temperature: Potential temperature for building a fixed hydrostatically-balanced reference state used in base-state subtraction. Can be a constant θ₀ or a function θ(z). Default: nothing (no base-state correction). When provided, an ExnerReferenceState is built during materialization.

struct ExplicitTimeStepping

Standard explicit time discretization for compressible dynamics.

All tendencies (including pressure gradient and acoustic modes) are computed together and time-stepped explicitly. This requires small time steps limited by the acoustic CFL condition (sound speed ~340 m/s).

Use SplitExplicitTimeDiscretization for more efficient time-stepping with larger Δt.

struct SplitExplicitTimeDiscretization{N, FT}

Split-explicit time discretization for compressible dynamics using the Exner pressure formulation following CM1 (Bryan 2002).

Uses acoustic substepping following Wicker and Skamarock (2002):

• Outer loop: WS-RK3 for slow tendencies (advection, Coriolis, diffusion)
• Inner loop: Forward-backward acoustic substeps for fast tendencies (pressure gradient)

The acoustic loop uses velocity (u, v, w) and Exner pressure perturbation (π') as prognostic variables, with a vertically implicit w-π' coupling and CM1-style divergence damping.

This allows using advective CFL time steps (~10-20 m/s) instead of acoustic CFL time steps (~340 m/s), typically enabling ~6x larger time steps.

Fields

• substeps: Number of acoustic substeps for the full time step (stage 3 of WS-RK3). For WS-RK3, earlier stages take fewer substeps ($Nτ = \mathrm{round}(β N)$), keeping $Δτ = Δt/N$ constant. Default: nothing (automatically computed from the acoustic CFL condition each time step)
• forward_weight: Off-centering parameter ω for the vertically implicit solver. ω > 0.5 damps vertical acoustic modes. Default: 0.6 (CM1 default)
• divergence_damping_coefficient: Forward-extrapolation filter coefficient $ϰ^{di}$ applied to the Exner pressure perturbation: $π̃' = π' + ϰ^{di} (π' - π'_{old})$. Default: 0.10 (CM1 default)
• acoustic_damping_coefficient: Klemp (2018) divergence damping $ϰ^{ac}$. Post-implicit-solve velocity correction: $u -= ϰ^{ac} c_p θ_v ∂Δπ'/∂x$. Provides constant damping per outer Δt regardless of substep count. Needed by WS-RK3 at large Δt. Default: 0.0

See also ExplicitTimeStepping.

acoustic_rk3_substep_loop!(
    model,
    substepper,
    Δt,
    β_stage,
    U⁰
)

Execute the acoustic substep loop for a Wicker-Skamarock RK3 stage using the Exner pressure formulation.

The acoustic substep size is constant: $Δτ = Δt / N$. Each stage takes $Nτ = \max(\mathrm{round}(β N), 1)$ substeps.

acoustic_substep_loop!(model, substepper, Δt, α_ssp, U⁰)

Execute the acoustic substep loop for an SSP RK3 stage. Delegates to the same Exner pressure acoustic loop used by WS-RK3.

prepare_acoustic_cache!(substepper, model)

Prepare the acoustic cache for an RK stage.

Computes stage-frozen coefficients for the Exner pressure acoustic loop:

1. Virtual potential temperature θᵥ (frozen during acoustic loop)
2. Pressure tendency coefficient S = c²/(cᵖ ρ₀ θᵥ²)
3. Exner pressure perturbation π' = (p/pˢᵗ)^(R/cᵖ) - π₀
4. Reference Exner function π₀ from the reference state

Following CM1's sound.F, the acoustic loop prognostics velocity and Exner pressure perturbation, with density diagnosed from the equation of state after the loop.

SubsidenceForcing(
    wˢ
) -> SubsidenceForcing{_A, Nothing, Nothing} where _A

Forcing that represents large-scale subsidence advecting horizontally-averaged fields downward:

\[F_{ρ ϕ} = - ρᵣ wˢ ∂_z \overline{ϕ}\]

where $wˢ$ is the subsidence_vertical_velocity, $ρᵣ$ is the reference density, and $\overline{ϕ}$ is the horizontal average of the field being forced.

Fields

• wˢ: Either a function of z specifying the subsidence velocity profile, or a Field containing the subsidence velocity.

The horizontal average is computed automatically during update_state!.

Example

using Breeze

grid = RectilinearGrid(size=(64, 64, 75), x=(0, 6400), y=(0, 6400), z=(0, 3000))

wˢ(z) = z < 1500 ? -0.0065 * z / 1500 : -0.0065 * (1 - (z - 1500) / 600)
subsidence = SubsidenceForcing(wˢ)
forcing = (; ρθ=subsidence, ρqᵗ=subsidence)

model = AtmosphereModel(grid; forcing)

model.forcing.ρθ

# output
SubsidenceForcing with wˢ: 1×1×76 Field{Nothing, Nothing, Face} reduced over dims = (1, 2) on RectilinearGrid on CPU
└── averaged_field: 1×1×75 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on RectilinearGrid on CPU

geostrophic_forcings(
    uᵍ,
    vᵍ
) -> NamedTuple{(:ρu, :ρv), <:Tuple{Breeze.Forcings.GeostrophicForcing{Oceananigans.Grids.XDirection, _A, Nothing} where _A, Breeze.Forcings.GeostrophicForcing{Oceananigans.Grids.YDirection, _A, Nothing} where _A}}

Create a pair of geostrophic forcings for the x- and y-momentum equations.

The Coriolis parameter is extracted from the model's coriolis during model construction.

Arguments

• uᵍ: Function of z specifying the x-component of the geostrophic velocity.
• vᵍ: Function of z specifying the y-component of the geostrophic velocity.

Returns a NamedTuple with ρu and ρv forcing entries that can be merged into the model forcing.

Example

using Breeze

uᵍ(z) = -10 + 0.001z
vᵍ(z) = 0.0

coriolis = FPlane(f=1e-4)
forcing = geostrophic_forcings(uᵍ, vᵍ)

# output
NamedTuple with 2 GeostrophicForcings:
├── ρu: GeostrophicForcing{XDirection}
│   └── geostrophic_momentum: vᵍ (generic function with 1 method)
└── ρv: GeostrophicForcing{YDirection}
    └── geostrophic_momentum: uᵍ (generic function with 1 method)

KinematicDriver

Module implementing kinematic dynamics for atmosphere models.

Kinematic dynamics prescribes the velocity field rather than solving for it, enabling isolated testing of microphysics, thermodynamics, and other physics without the complexity of solving the momentum equations.

This is analogous to the kin1d driver in P3-microphysics.

Wrapper indicating that density is fixed (not prognostic).

struct PrescribedDynamics{Div, D, P, FT}

Dynamics for kinematic atmosphere models where velocity is prescribed. The type parameter Div indicates whether divergence correction is applied.

PrescribedDynamics(
    density;
    pressure,
    surface_pressure,
    standard_pressure,
    divergence_correction
) -> PrescribedDynamics{_A, D} where {_A, D<:PrescribedDensity}

Construct PrescribedDynamics from a density field or PrescribedDensity. If pressure=nothing, hydrostatic pressure is computed during materialization.

PrescribedDynamics(
    reference_state::ReferenceState;
    divergence_correction
) -> PrescribedDynamics{_A, D} where {_A, D<:PrescribedDensity}

Construct PrescribedDynamics from a ReferenceState. Wraps density in PrescribedDensity (fixed in time).

If divergence_correction=true, scalar tendencies include +c∇·(ρU) to account for the non-divergent velocity field.

Example

using Oceananigans
using Breeze

grid = RectilinearGrid(size=(4, 4, 8), extent=(1000, 1000, 2000))
reference_state = ReferenceState(grid, ThermodynamicConstants())
dynamics = PrescribedDynamics(reference_state)

# output
PrescribedDynamics
├── density: PrescribedDensity
├── pressure: 1×1×8 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on RectilinearGrid on CPU
├── surface_pressure: 101325.0
└── standard_pressure: 100000.0

BulkMicrophysics(
;
    ...
) -> BulkMicrophysics{N, Nothing, Nothing} where N<:(SaturationAdjustment{E} where E<:MixedPhaseEquilibrium)
BulkMicrophysics(
    FT::DataType;
    categories,
    cloud_formation,
    precipitation_boundary_condition
) -> BulkMicrophysics{N, Nothing, Nothing} where N<:(SaturationAdjustment{E} where E<:MixedPhaseEquilibrium)

Return a BulkMicrophysics microphysics scheme.

Keyword arguments

• categories: Microphysical categories (e.g., cloud liquid, cloud ice, rain, snow) or nothing for non-precipitating
• cloud_formation: Cloud formation scheme (default: SaturationAdjustment)
• precipitation_boundary_condition: Bottom boundary condition for precipitation sedimentation.
  • nothing (default): Precipitation passes through the bottom
  • ImpenetrableBoundaryCondition(): Precipitation collects at the bottom

struct BulkMicrophysics{N, C, B}

Bulk microphysics scheme with cloud formation and precipitation categories.

Fields

• cloud_formation: Cloud formation scheme (saturation adjustment or non-equilibrium)
• categories: Precipitation categories (e.g., rain, snow) or nothing
• precipitation_boundary_condition: Bottom boundary condition for precipitation sedimentation.
  • nothing (default): Precipitation passes through the bottom (open boundary)
  • ImpenetrableBoundaryCondition(): Precipitation collects at the bottom (zero terminal velocity at surface)

Return a condensate formation model that applies a constant phase-change rate.

This type is intended to be usable for both liquid (condensation/evaporation) and ice (deposition/sublimation).

DCMIP2016KesslerMicrophysics(
;
    ...
) -> Breeze.Microphysics.DCMIP2016KesslerMicrophysics
DCMIP2016KesslerMicrophysics(
    FT;
    dcmip_temperature_scale,
    terminal_velocity_coefficient,
    density_scale,
    terminal_velocity_exponent,
    autoconversion_rate,
    autoconversion_threshold,
    accretion_rate,
    accretion_exponent,
    evaporation_ventilation_coefficient_1,
    evaporation_ventilation_coefficient_2,
    evaporation_ventilation_exponent_1,
    evaporation_ventilation_exponent_2,
    diffusivity_coefficient,
    thermal_conductivity_coefficient,
    substep_cfl
) -> Breeze.Microphysics.DCMIP2016KesslerMicrophysics

Construct a DCMIP2016 implementation of the Kessler (1969) warm-rain bulk microphysics scheme.

This implementation follows the DCMIP2016 test case specification, which is based on Klemp and Wilhelmson (1978).

Positional Arguments

• FT: Floating-point type for all parameters (default: Oceananigans.defaults.FloatType).

References

• Zarzycki, C. M., et al. (2019). DCMIP2016: the splitting supercell test case. Geoscientific Model Development, 12, 879–892.
• Kessler, E. (1969). On the Distribution and Continuity of Water Substance in Atmospheric Circulations. Meteorological Monographs, 10(32).
• Klemp, J. B., & Wilhelmson, R. B. (1978). The Simulation of Three-Dimensional Convective Storm Dynamics. Journal of the Atmospheric Sciences, 35(6), 1070-1096.
• DCMIP2016 Fortran implementation (kessler.f90 in DOI: 10.5281/zenodo.1298671)

Moisture Categories

This scheme represents moisture in three categories:

• Water vapor mixing ratio (rᵛ)
• Cloud water mixing ratio (rᶜˡ)
• Rain water mixing ratio (rʳ)

Breeze tracks moisture using mass fractions (q), whereas the Kessler scheme uses mixing ratios (r). Conversions between these representations are performed internally. In Breeze, water vapor is not a prognostic variable; instead, it is diagnosed from the total specific moisture qᵗ and the liquid condensates.

Physical Processes

1. Autoconversion: Cloud water converts to rain water when the cloud water mixing ratio exceeds a threshold.
2. Accretion: Rain water collects cloud water as it falls.
3. Saturation Adjustment: Water vapor condenses to cloud water or cloud water evaporates to maintain saturation.
4. Rain Evaporation: Rain water evaporates into subsaturated air.
5. Rain Sedimentation: Rain water falls gravitationally.

Implementation Details

• The microphysics update is applied via a GPU-compatible kernel launched from microphysics_model_update!.
• Rain sedimentation uses subcycling to satisfy CFL constraints, following the Fortran implementation.
• All microphysical updates are applied directly to the state variables in the kernel.

Keyword Arguments

Saturation (Tetens/Clausius-Clapeyron formula)

• dcmip_temperature_scale (T_DCMIP2016): A parameter of uncertain provenance that appears in the DCMIP2016 implementation of the Kessler scheme (line 105 of kessler.f90 in DOI: 10.5281/zenodo.1298671)

The "saturation adjustment coefficient" f₅ is then computed as

\[f₅ = a × T_DCMIP2016 × ℒˡᵣ / cᵖᵈ\]

where a is the liquid_coefficient for Tetens' saturation vapor pressure formula, ℒˡᵣ is the latent heat of vaporization of liquid water, and cᵖᵈ is the heat capacity of dry air.

Rain Terminal Velocity (Klemp & Wilhelmson 1978, eq. 2.15)

Terminal velocity: 𝕎ʳ = a𝕎 × (ρ × rʳ × Cᵨ)^β𝕎 × √(ρ₀/ρ)

• terminal_velocity_coefficient (a𝕎): Terminal velocity coefficient in m/s (default: 36.34)
• density_scale (Cᵨ): Density scale factor for unit conversion (default: 0.001)
• terminal_velocity_exponent (β𝕎): Terminal velocity exponent (default: 0.1364)
• ρ: Density
• ρ₀: Reference density at z=0

Autoconversion

• autoconversion_rate (k₁): Autoconversion rate coefficient in s⁻¹ (default: 0.001)
• autoconversion_threshold (rᶜˡ★): Critical cloud water mixing ratio threshold in kg/kg (default: 0.001)

Accretion

• accretion_rate (k₂): Accretion rate coefficient in s⁻¹ (default: 2.2)
• accretion_exponent (βᵃᶜᶜ): Accretion exponent for rain mixing ratio (default: 0.875)

Rain Evaporation (Klemp & Wilhelmson 1978, eq. 2.14)

Ventilation: (Cᵉᵛ₁ + Cᵉᵛ₂ × (ρ rʳ)^βᵉᵛ₁) × (ρ rʳ)^βᵉᵛ₂

• evaporation_ventilation_coefficient_1 (Cᵉᵛ₁): Evaporation ventilation coefficient 1 (default: 1.6)
• evaporation_ventilation_coefficient_2 (Cᵉᵛ₂): Evaporation ventilation coefficient 2 (default: 124.9)
• evaporation_ventilation_exponent_1 (βᵉᵛ₁): Evaporation ventilation exponent 1 (default: 0.2046)
• evaporation_ventilation_exponent_2 (βᵉᵛ₂): Evaporation ventilation exponent 2 (default: 0.525)
• diffusivity_coefficient (Cᵈⁱᶠᶠ): Diffusivity-related denominator coefficient (default: 2.55e8)
• thermal_conductivity_coefficient (Cᵗʰᵉʳᵐ): Thermal conductivity-related denominator coefficient (default: 5.4e5)

Numerical

• substep_cfl: CFL safety factor for sedimentation subcycling (default: 0.8)

struct DCMIP2016KesslerMicrophysics{FT}

DCMIP2016 implementation of the Kessler (1969) warm-rain bulk microphysics scheme. See the constructor DCMIP2016KesslerMicrophysics for full documentation.

NonEquilibriumCloudFormation(liquid, ice=nothing)

A cloud formation scheme where cloud liquid and ice are prognostic variables that evolve via condensation/evaporation and deposition/sublimation tendencies, rather than being diagnosed instantaneously via saturation adjustment.

The condensation/evaporation and deposition/sublimation tendencies are commonly modeled as relaxation toward saturation with timescale τ_relax, including a latent-heat (psychrometric/thermal) correction factor; see Morrison and Grabowski (2008), Appendix Eq. (A3), and standard cloud microphysics texts such as Pruppacher and Klett (2010) or Rogers and Yau (1989).

For some bulk schemes (e.g. the CloudMicrophysics 1M extension), liquid and ice may be set to nothing and used purely as phase indicators (warm-phase vs mixed-phase), with any relaxation timescales sourced from the scheme's precipitation/category parameters instead.

Fields

• liquid: Parameters for cloud liquid (contains relaxation timescale τ_relax)
• ice: Parameters for cloud ice (contains relaxation timescale τ_relax), or nothing for warm-phase only

References

• Morrison, H. and Grabowski, W. W. (2008). A novel approach for representing ice microphysics in models: Description and tests using a kinematic framework. J. Atmos. Sci., 65, 1528–1548. https://doi.org/10.1175/2007JAS2491.1
• Pruppacher, H. R. and Klett, J. D. (2010). Microphysics of Clouds and Precipitation (2nd ed.).
• Rogers, R. R. and Yau, M. K. (1989). A Short Course in Cloud Physics (3rd ed.).

SaturationAdjustment(
;
    ...
) -> SaturationAdjustment{E} where E<:MixedPhaseEquilibrium
SaturationAdjustment(
    FT::DataType;
    tolerance,
    maxiter,
    equilibrium
) -> SaturationAdjustment{E} where E<:MixedPhaseEquilibrium

Return SaturationAdjustment microphysics representing an instantaneous adjustment to equilibrium between condensates and water vapor, computed by a solver with tolerance and maxiter.

The options for equilibrium are:

• WarmPhaseEquilibrium() representing an equilibrium between water vapor and liquid water.

• MixedPhaseEquilibrium() representing a temperature-dependent equilibrium between water vapor, possibly supercooled liquid water, and ice. The equilibrium state is modeled as a linear variation of the equilibrium liquid fraction with temperature, between the freezing temperature (e.g. 273.15 K) below which liquid water is supercooled, and the temperature of homogeneous ice nucleation temperature (e.g. 233.15 K) at which the supercooled liquid fraction vanishes.

RelativeHumidity(
    model
) -> KernelFunctionOperation{Center, Center, Center, _A, _B, K, Tuple{}} where {_A, _B, K<:Breeze.Microphysics.RelativeHumidityKernelFunction}

Return a KernelFunctionOperation representing the relative humidity $ℋ$, defined as the ratio of vapor pressure to saturation vapor pressure:

\[ℋ = \frac{pᵛ}{pᵛ⁺}\]

where $pᵛ$ is the vapor pressure (partial pressure of water vapor) computed from the ideal gas law

\[pᵛ = ρ qᵛ Rᵛ T\]

and $pᵛ⁺$ is the saturation vapor pressure.

For unsaturated conditions, $ℋ < 1$. For saturated conditions with saturation adjustment microphysics, $ℋ = 1$ (or very close to it due to numerical precision).

Examples

using Breeze
grid = RectilinearGrid(size=(1, 1, 128), extent=(1e3, 1e3, 1e3))
microphysics = SaturationAdjustment()
model = AtmosphereModel(grid; microphysics)
set!(model, θ=300, qᵗ=0.005)  # subsaturated
ℋ = RelativeHumidity(model)

# output
KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── kernel_function: RelativeHumidityKernelFunction
└── arguments: ()

As with other diagnostics, RelativeHumidity may be wrapped in Field to store the result:

ℋ_field = RelativeHumidity(model) |> Field

# output
1×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×134 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:131) with eltype Float64 with indices 0:2×0:2×-2:131
    └── max=0.214947, min=0.136946, mean=0.172492

We also provide a convenience constructor for the Field:

ℋ_field = RelativeHumidityField(model)

# output
1×1×128 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 3×3×134 OffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:131) with eltype Float64 with indices 0:2×0:2×-2:131
    └── max=0.214947, min=0.136946, mean=0.172492

adjust_thermodynamic_state(
    𝒰₀::Breeze.Thermodynamics.AbstractThermodynamicState,
    microphysics::SaturationAdjustment,
    constants
) -> Breeze.Thermodynamics.AbstractThermodynamicState

Return the saturation-adjusted thermodynamic state using a secant iteration.

compute_temperature(
    𝒰₀,
    adjustment::SaturationAdjustment,
    constants
) -> Any

Perform saturation adjustment and return the temperature associated with the adjusted state.

kessler_terminal_velocity(rʳ, ρ, ρ₁, microphysics) -> Any

Compute rain terminal velocity (m/s) following Klemp and Wilhelmson (1978) eq. 2.15.

The terminal velocity is computed as:

\[𝕎ʳ = a^𝕎 × (ρ × rʳ × Cᵨ)^{β^𝕎} × \sqrt{ρ₀/ρ}\]

where a^𝕎 is terminal_velocity_coefficient, Cᵨ is density_scale, and β^𝕎 is terminal_velocity_exponent.

MoistAirBuoyancy(
    grid;
    surface_pressure,
    reference_potential_temperature,
    standard_pressure,
    thermodynamic_constants
) -> MoistAirBuoyancy{RS, AT} where {RS<:(ReferenceState{_A, P, D, T, QV, QL, QI} where {_A, P<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), D<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), T<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), QV<:(Oceananigans.Fields.ZeroField{T, 3} where T), QL<:(Oceananigans.Fields.ZeroField{T, 3} where T), QI<:(Oceananigans.Fields.ZeroField{T, 3} where T)}), AT<:ThermodynamicConstants}

Return a MoistAirBuoyancy formulation that can be provided as input to an Oceananigans.NonhydrostaticModel.

Example

using Breeze, Oceananigans

grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 3e3))
buoyancy = MoistAirBuoyancy(grid)

# output
MoistAirBuoyancy:
├── reference_state: ReferenceState{Float64}(p₀=101325.0, θ₀=288.0, pˢᵗ=100000.0)
└── thermodynamic_constants: ThermodynamicConstants{Float64}

To build a model with MoistAirBuoyancy, we include potential temperature and total specific humidity tracers θ and qᵗ to the model.

model = NonhydrostaticModel(grid; buoyancy, tracers = (:θ, :qᵗ))

# output
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (θ, qᵗ)
├── closure: Nothing
├── buoyancy: MoistAirBuoyancy with ĝ = NegativeZDirection()
└── coriolis: Nothing

ParcelDynamics(
;
    ...
) -> ParcelDynamics{Nothing, Nothing, Nothing, Nothing}
ParcelDynamics(
    FT::DataType;
    surface_pressure,
    standard_pressure
) -> ParcelDynamics{Nothing, Nothing, Nothing, Nothing}

Construct ParcelDynamics with default (uninitialized) state.

The environmental profiles and parcel state are set using set! after constructing the AtmosphereModel.

struct ParcelDynamics{S, TS, D, P, FT}

Lagrangian parcel dynamics for AtmosphereModel.

Fields

• state: parcel state (position, thermodynamics, microphysics)
• timestepper: SSP RK3 timestepper with tendencies
• density: environmental density field [kg/m³]
• pressure: environmental pressure field [Pa]
• surface_pressure: surface pressure [Pa]
• standard_pressure: standard pressure for potential temperature [Pa]

mutable struct ParcelInitialState{FT, MP}

Storage for the initial parcel prognostic state at the beginning of a time step. Used by SSP RK3 to combine the initial state with intermediate states.

Fields

• x, y, z: initial position [m]
• qᵗ: initial specific total moisture [kg/kg]
• ℰ: initial specific static energy [J/kg] or potential temperature [K]
• μ: initial microphysics prognostics (density-weighted)

ParcelModel

Type alias for AtmosphereModel{<:ParcelDynamics}.

A ParcelModel represents a Lagrangian adiabatic parcel that rises through a prescribed environmental atmosphere. The parcel is characterized by its position (x, y, z), thermodynamic state, and moisture content. The environmental profiles (temperature, pressure, density, velocities) are defined on a 1D vertical grid.

The parcel's motion is determined by interpolating environmental velocities to the parcel position, and its thermodynamic evolution follows adiabatic processes with optional microphysical interactions.

See also ParcelDynamics, AtmosphereModel.

mutable struct ParcelState{FT, TH, MP}

• x::Any

• y::Any

• z::Any

• ρ::Any

• qᵗ::Any

• ρqᵗ::Any

• ℰ::Any

• ρℰ::Any

• 𝒰::Any

• μ::Any

State of a Lagrangian air parcel with position, thermodynamic state, and microphysics.

The parcel model evolves specific quantities (qᵗ, ℰ) directly for exact conservation. Density-weighted forms (ρqᵗ, ρℰ) are also stored for consistency with the microphysics interface.

mutable struct ParcelTendencies{FT, GM}

Tendencies (time derivatives) for parcel prognostic variables.

Fields

• Gx, Gy, Gz: position tendencies [m/s]
• Ge: specific energy tendency [J/kg/s]
• Gqᵗ: specific moisture tendency [kg/kg/s]
• Gμ: microphysics prognostic tendencies (density-weighted)

struct ParcelTimestepper{GT, U0, FT}

SSP RK3 time-stepper for ParcelModel.

Stores tendencies, the initial state at the beginning of a time step, and the SSP RK3 stage coefficients.

Fields

• G: tendencies for prognostic variables
• U⁰: initial state storage (position, moisture, thermodynamics, microphysics)
• α¹, α², α³: SSP RK3 stage coefficients (1, 1/4, 2/3)

ParcelTimestepper(
    state::ParcelState{FT},
    Gμ
) -> Breeze.ParcelModels.ParcelTimestepper{GT, U0} where {GT<:Breeze.ParcelModels.ParcelTendencies, U0<:Breeze.ParcelModels.ParcelInitialState}

Construct a ParcelTimestepper for SSP RK3 time-stepping.

Adjust the thermodynamic state for adiabatic ascent/descent to a new height. Conserves the thermodynamic variable (static energy or potential temperature).

compute_parcel_tendencies!(
    model::AtmosphereModel{<:ParcelDynamics}
)

Compute tendencies for the parcel prognostic variables.

Position tendencies are interpolated from environmental velocity fields. Thermodynamic and moisture tendencies come from microphysical sources/sinks.

The parcel model evolves specific quantities (e, qᵗ) directly, not density-weighted quantities. For adiabatic ascent with no microphysics, specific static energy and moisture are exactly conserved (de/dt = dqᵗ/dt = 0). This is simpler and more accurate than stepping density-weighted quantities.

materialize_parcel_microphysics_prognostics(
    FT,
    microphysics
) -> Union{Nothing, NamedTuple}

Create the parcel microphysics prognostic variables for the given microphysics scheme.

Returns nothing for microphysics schemes without explicit prognostic variables (e.g., Nothing, SaturationAdjustment), or a NamedTuple containing the prognostic density-weighted scalars for schemes with prognostic microphysics.

The prognostic variables use the same ρ-weighted names as the grid-based model (e.g., :ρqᶜˡ, :ρqʳ) from prognostic_field_names(microphysics).

ssp_rk3_parcel_substep!(
    model::AtmosphereModel{<:ParcelDynamics},
    U⁰::Breeze.ParcelModels.ParcelInitialState,
    Δt,
    α
)

Apply an SSP RK3 substep with coefficient α:

\[u^{(m)} = (1 - α) u^{(0)} + α (u^{(m-1)} + Δt G^{(m-1)})\]

where u^{(0)} is the initial state, u^{(m-1)} is the current state, and G^{(m-1)} is the tendency at the current state.

The parcel model steps specific quantities (e, qᵗ) directly for exact conservation. For adiabatic ascent with no microphysics sources, de/dt = dqᵗ/dt = 0, so these quantities remain exactly constant throughout the simulation.

step_parcel_state!(
    model::AtmosphereModel{<:ParcelDynamics},
    Δt
)

Step the parcel state forward using Forward Euler: x^(n+1) = x^n + Δt * G^n.

Computes tendencies at the current state, then advances all prognostic variables. After updating position, the thermodynamic state is adjusted for the new height (adiabatic adjustment) and environmental conditions are updated from the profiles.

store_initial_parcel_state!(
    U⁰::Breeze.ParcelModels.ParcelInitialState,
    state::ParcelState
)

Copy current prognostic state values to the initial state storage.

PotentialTemperatureFormulations

Submodule defining the liquid-ice potential temperature thermodynamic formulation for atmosphere models.

LiquidIcePotentialTemperatureFormulation uses liquid-ice potential temperature density ρθ as the prognostic thermodynamic variable.

LiquidIcePotentialTemperatureFormulation uses liquid-ice potential temperature density ρθ as the prognostic thermodynamic variable.

Liquid-ice potential temperature is a conserved quantity in moist adiabatic processes and is defined as:

\[θˡⁱ = T \left( \frac{p^{st}}{p} \right)^{Rᵐ/cᵖᵐ} \exp\left( \frac{ℒˡᵣ qˡ + ℒⁱᵣ qⁱ}{cᵖᵐ T} \right)\]

StaticEnergyFormulations

Submodule defining the static energy thermodynamic formulation for atmosphere models.

StaticEnergyFormulation uses moist static energy density ρe as the prognostic thermodynamic variable. Moist static energy is a conserved quantity in adiabatic, frictionless flow that combines sensible heat, gravitational potential energy, and latent heat.

StaticEnergyFormulation uses moist static energy density ρe as the prognostic thermodynamic variable.

Moist static energy is a conserved quantity in adiabatic, frictionless flow that combines sensible heat, gravitational potential energy, and latent heat:

\[e = cᵖᵐ T + g z - ℒˡᵣ qˡ - ℒⁱᵣ qⁱ\]

The energy density equation includes a buoyancy flux term following Pauluis (2008).

struct ClausiusClapeyron

A saturation vapor pressure formulation based on the Clausius-Clapeyron relation.

The Clausius-Clapeyron equation describes how saturation vapor pressure varies with temperature based on thermodynamic principles. This formulation uses thermodynamic constants (latent heats, heat capacities, triple point values) to compute saturation vapor pressure analytically.

See saturation_vapor_pressure for the implementation details.

ClausiusClapeyronThermodynamicConstants{FT, C, I}

Type alias for ThermodynamicConstants using the Clausius-Clapeyron formulation for saturation vapor pressure calculations.

CondensedPhase(; ...)
CondensedPhase(FT; reference_latent_heat, heat_capacity)

Return CondensedPhase with specified parameters converted to FT.

Two examples of CondensedPhase are liquid and ice. When matter is converted from vapor to liquid, water molecules in the gas phase cluster together and slow down to form liquid with heat_capacity, The lost of molecular kinetic energy is called the reference_latent_heat.

Likewise, during deposition, water molecules in the gas phase cluster into ice crystals.

Arguments

• FT: Float type to use (defaults to Oceananigans.defaults.FloatType)
• reference_latent_heat: Difference between the internal energy of the gaseous phase at the energy_reference_temperature.
• heat_capacity: Heat capacity of the phase of matter.

ExnerReferenceState(
    grid;
    ...
) -> ExnerReferenceState{_A, FP, FD, FE} where {_A, FP<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), FD<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), FE<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B})}
ExnerReferenceState(
    grid,
    constants;
    surface_pressure,
    potential_temperature,
    standard_pressure
) -> ExnerReferenceState{_A, FP, FD, FE} where {_A, FP<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), FD<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), FE<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B})}

Construct an ExnerReferenceState by discrete Exner integration on grid.

The Exner function π₀ is built by integrating upward from the surface:

\[π₀[k] = π₀[k-1] - \frac{g \, Δz}{cᵖᵈ \, θᵣ^{face}[k]}\]

where $θᵣ^{face}$ is the face-averaged reference potential temperature. This ensures the discrete Exner hydrostatic balance is exact.

Pressure and density are then derived from π₀:

• $p₀ = pˢᵗ \, π₀^{cᵖᵈ/Rᵈ}$
• $T₀ = θᵣ \, π₀$
• $ρ₀ = p₀ / (Rᵈ \, T₀)$

Arguments

• grid: The grid
• constants: Thermodynamic constants (default: ThermodynamicConstants(eltype(grid)))

Keyword Arguments

• surface_pressure: Pressure at z=0 (default: 101325 Pa)
• potential_temperature: Constant value or function θᵣ(z) for reference potential temperature (default: 288 K)
• standard_pressure: pˢᵗ for potential temperature definition (default: 1e5 Pa)

ExnerReferenceState

A dry reference state built in Exner coordinates, ensuring that the discrete Exner hydrostatic balance

\[cᵖᵈ θᵣ^{face} \frac{π₀[k] - π₀[k-1]}{Δz} = -g\]

holds EXACTLY at every interior z-face. This is essential for the Exner pressure acoustic substepping formulation, where the vertical pressure gradient is computed as $cᵖᵈ θᵥ ∂π'/∂z$ and the hydrostatic part must cancel to machine precision.

Unlike ReferenceState which builds pressure first and derives Exner, this type builds the Exner function π₀ first by discrete integration and then derives pressure and density from it. This matches CM1's approach where pi0 is the fundamental reference variable.

Fields

• surface_pressure: Reference pressure at z=0 (Pa)
• surface_potential_temperature: Reference potential temperature at z=0 (K)
• standard_pressure: pˢᵗ for potential temperature definition (Pa)
• pressure: Reference pressure field $p₀ = pˢᵗ π₀^{cᵖᵈ/Rᵈ}$ (derived from π₀)
• density: Reference density field $ρ₀ = p₀/(Rᵈ T₀)$ (derived from π₀ and θᵣ)
• exner_function: Reference Exner function π₀ (built by discrete integration)

struct IdealGas{FT}

A struct representing an ideal gas with molar mass and specific heat capacity.

Fields

• molar_mass: Molar mass of the gas in kg/mol
• heat_capacity: Specific heat capacity at constant pressure in J/(kg·K)

Examples

using Breeze
dry_air = IdealGas(molar_mass=0.02897, heat_capacity=1005)

# output
IdealGas{Float64}(molar_mass=0.02897, heat_capacity=1005.0)

MixedPhaseEquilibrium(; freezing_temperature=273.15, homogeneous_ice_nucleation_temperature=233.15)

Represents a mixed-phase equilibrium where both liquid and ice condensates are considered. The liquid fraction varies linearly with temperature between the freezing temperature and the homogeneous ice nucleation temperature.

struct MoistureMassFractions{FT}

A struct representing the moisture mass fractions of a moist air parcel.

Fields

• vapor: the mass fraction of vapor
• liquid: the mass fraction of liquid
• ice: the mass fraction of ice

MoistureMassFractions(
    r::Breeze.Thermodynamics.MoistureMixingRatio
) -> Breeze.Thermodynamics.MoistureMassFractions

Convert MoistureMixingRatio to MoistureMassFractions.

Mass fractions are defined as mass of constituent per total mass:

\[q = r / (1 + rᵗ)\]

where rᵗ is the total mixing ratio.

MoistureMixingRatio(
    q::Breeze.Thermodynamics.MoistureMassFractions
) -> Breeze.Thermodynamics.MoistureMixingRatio

Convert MoistureMassFractions to MoistureMixingRatio.

Mixing ratios are defined as mass of constituent per mass of dry air:

\[r = q / (1 - qᵗ) = q / qᵈ\]

where qᵗ is the total specific moisture and qᵈ = 1 - qᵗ is the dry air mass fraction.

struct PlanarMixedPhaseSurface{FT}

Return PlanarMixedPhaseSurface for computing the saturation vapor pressure over a surface composed of a mixture of liquid and ice, with a given liquid_fraction.

ReferenceState(
    grid;
    ...
) -> ReferenceState{_A, P, D, T, QV, QL, QI} where {_A, P<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), D<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), T<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), QV<:(Oceananigans.Fields.ZeroField{T, 3} where T), QL<:(Oceananigans.Fields.ZeroField{T, 3} where T), QI<:(Oceananigans.Fields.ZeroField{T, 3} where T)}
ReferenceState(
    grid,
    constants;
    surface_pressure,
    potential_temperature,
    standard_pressure,
    discrete_hydrostatic_balance,
    vapor_mass_fraction,
    liquid_mass_fraction,
    ice_mass_fraction
) -> ReferenceState{_A, P, D, T, QV, QL, QI} where {_A, P<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), D<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), T<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), QV<:(Oceananigans.Fields.ZeroField{T, 3} where T), QL<:(Oceananigans.Fields.ZeroField{T, 3} where T), QI<:(Oceananigans.Fields.ZeroField{T, 3} where T)}

Return a ReferenceState on grid, with ThermodynamicConstants constants that includes the hydrostatic reference pressure and reference density.

The reference state is initialized with a dry adiabatic temperature profile and the given moisture profiles (zero by default). The pressure and density are then computed by hydrostatic integration using the mixture gas constant Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ and the ideal gas law ρ = p / (Rᵐ T).

Arguments

• grid: The grid.
• constants :: ThermodynamicConstants: By default, ThermodynamicConstants(eltype(grid)).

Keyword arguments

• surface_pressure: By default, 101325.

• potential_temperature: A constant value (default 288) or a function θ(z) giving the potential temperature profile. When a constant is provided, closed-form adiabatic hydrostatic profiles are used. When a function is provided, the hydrostatic profiles are computed by numerical integration of dp/dz = -gρ.

• standard_pressure: Reference pressure for potential temperature (pˢᵗ). By default, 1e5.

• discrete_hydrostatic_balance: If true, recompute the reference pressure from the reference density using discrete integration, so that ∂z(p_ref) + g * ℑz(ρ_ref) = 0 exactly at the discrete level. By default, false.

  Discrete vs continuous hydrostatic balance

  With discrete balance, reference subtraction becomes a no-op (the subtracted terms cancel to machine precision). For split-explicit compressible dynamics, continuous balance (default) is preferred: both the actual and reference states share similar $O(Δz^2)$ truncation errors that cancel in the perturbation PG, leaving only the tiny truncation error of the physical perturbation $∂(p - p_{ref})/∂z$.

• vapor_mass_fraction: Initial qᵛ profile. Can be a Number, Function(z), or Field. Default: nothing (ZeroField).

• liquid_mass_fraction: Initial qˡ profile. Default: nothing (ZeroField).

• ice_mass_fraction: Initial qⁱ profile. Default: nothing (ZeroField).

Pass =0 to allocate an actual Field initialized to zero — required for later use with compute_reference_state! or set_to_mean!.

TetensFormula(; ...) -> Breeze.Thermodynamics.TetensFormula
TetensFormula(
    FT;
    reference_saturation_vapor_pressure,
    reference_temperature,
    liquid_coefficient,
    liquid_temperature_offset,
    ice_coefficient,
    ice_temperature_offset
) -> Breeze.Thermodynamics.TetensFormula

Construct a TetensFormula saturation vapor pressure formulation. Tetens's (1930) formula is an empirical relationship for the saturation vapor pressure,

\[pᵛ⁺(T) = pᵛ⁺ᵣ \exp \left( a \frac{T - Tᵣ}{T - δT} \right) ,\]

where $pᵛ⁺ᵣ$ is reference_saturation_vapor_pressure, $Tᵣ$ is reference_temperature, $a$ is an empirical coefficient, and $δT$ is a temperature offset.

See also the wikipedia article on "Tetens equation". Different coefficients are used for liquid water and ice surfaces. Default values for the liquid formula are from Monteith and Unsworth (2014), and default values for the ice formula are from Murray (1967):

Liquid water (T > 0°C):

• liquid_coefficient: 17.27
• liquid_temperature_offset: 35.85 K (corresponding to 237.3 K offset from 0°C)

Ice (T < 0°C):

• ice_coefficient: 21.875
• ice_temperature_offset: 7.65 K (corresponding to 265.5 K offset from 0°C)

References

• Tetens, O. (1930). Über einige meteorologische Begriffe. Zeitschrift für Geophysik 6, 297–309.
• Monteith, J. L. and Unsworth, M. H. (2014). Principles of Environmental Physics. 4th Edition (Academic Press).
• Murray, F. W. (1967). On the computation of saturation vapor pressure. Journal of Applied Meteorology 6, 203–204.
• Wikipedia: Tetens equation; https://en.wikipedia.org/wiki/Tetens_equation

Example

julia> using Breeze.Thermodynamics

julia> tf = TetensFormula()
TetensFormula{Float64}(pᵣ=610.0, Tᵣ=273.15, aˡ=17.27, δTˡ=35.85, aⁱ=21.875, δTⁱ=7.65)

TetensFormulaThermodynamicConstants{FT, C, I}

Type alias for ThermodynamicConstants using the Tetens formula for saturation vapor pressure calculations.

ThermodynamicConstants(; ...) -> ThermodynamicConstants
ThermodynamicConstants(
    FT;
    molar_gas_constant,
    gravitational_acceleration,
    energy_reference_temperature,
    triple_point_temperature,
    triple_point_pressure,
    dry_air_molar_mass,
    dry_air_heat_capacity,
    vapor_molar_mass,
    vapor_heat_capacity,
    liquid,
    ice,
    saturation_vapor_pressure
) -> ThermodynamicConstants

Return ThermodynamicConstants with parameters that represent gaseous mixture of dry "air" and vapor, as well as condensed liquid and ice phases. The triple_point_temperature and triple_point_pressure may be combined with internal energy parameters for condensed phases to compute the vapor pressure at the boundary between vapor and a homogeneous sample of the condensed phase. The gravitational_acceleration parameter is included to compute ReferenceState quantities associated with hydrostatic balance.

WarmPhaseEquilibrium()

Represents a warm-phase equilibrium where only liquid water condensate is considered. The equilibrated surface is always a planar liquid surface.

adiabatic_hydrostatic_density(
    z,
    p₀,
    θ₀,
    pˢᵗ,
    constants
) -> Any

Compute the reference density at height z that associated with the reference pressure p₀, potential temperature θ₀, and standard pressure pˢᵗ. The reference density is defined as the density of dry air at the reference pressure and temperature.

adiabatic_hydrostatic_pressure(z, p₀, θ₀, constants) -> Any

Compute the reference pressure at height z that associated with the reference pressure p₀ and potential temperature θ₀. The reference pressure is defined as the pressure of dry air at the reference pressure and temperature.

adjustment_saturation_specific_humidity(
    T,
    pᵣ,
    qᵗ,
    constants,
    surface
) -> Any

Compute the saturation specific humidity $qᵛ⁺$ for use in saturation adjustment, assuming saturated conditions where condensate is present.

This function always uses the saturated formula (equation 37 in paper by Pressel et al. 2015):

\[qᵛ⁺ = ϵᵈᵛ (1 - qᵗ) \frac{pᵛ⁺}{pᵣ - pᵛ⁺}\]

where $ϵᵈᵛ = Rᵈ / Rᵛ ≈ 0.622$.

Unlike equilibrium_saturation_specific_humidity, this function does not check whether the air is actually saturated. It is intended for use within the saturation adjustment iteration where we assume saturated conditions throughout.

adjustment_saturation_specific_humidity(
    T,
    pᵣ,
    qᵗ,
    constants,
    equilibrium::Breeze.Thermodynamics.AbstractPhaseEquilibrium
) -> Any

Compute the adjustment saturation specific humidity using a phase equilibrium model to determine the condensation surface based on temperature T.

compute_hydrostatic_reference!(
    ref::ReferenceState,
    constants
)

Compute the hydrostatic reference pressure and density profiles from the temperature and moisture mass fraction profiles stored in ref.

The integration uses the mixture gas constant Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ (where qᵈ = 1 - qᵛ - qˡ - qⁱ) and the ideal gas law ρ = p / (Rᵐ T).

compute_reference_state!(
    ref::ReferenceState,
    T̄,
    q̄ᵗ,
    constants
)

Convenience method that assumes all moisture is vapor (no condensate in the reference state). Equivalent to compute_reference_state!(reference_state, T̄, q̄ᵗ, 0, 0, constants).

compute_reference_state!(
    ref::ReferenceState,
    T̄,
    q̄ᵛ,
    q̄ˡ,
    q̄ⁱ,
    constants
)

Recompute the reference pressure and density profiles by setting the reference temperature to T̄ and moisture mass fractions to q̄ᵛ, q̄ˡ, q̄ⁱ, then integrating the hydrostatic equation using the mixture gas constant Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ and ideal gas law ρ = p / (Rᵐ T).

T̄, q̄ᵛ, q̄ˡ, q̄ⁱ can be Numbers, Function(z)s, or Fields.

This function is useful for:

• Initialization: setting the reference state to match a non-constant-θ initial condition
• Runtime: calling from a callback to keep the reference state close to the evolving mean state

dewpoint_temperature(
    pᵛ,
    T,
    constants,
    surface;
    tolerance,
    maxiter
) -> Any

Compute the dewpoint temperature $T⁺$ given the vapor pressure pᵛ, actual temperature T, thermodynamic constants, and condensation surface.

The dewpoint temperature is defined as the temperature at which the saturation vapor pressure equals the actual vapor pressure:

\[pᵛ⁺(T⁺) = pᵛ\]

This implicit equation is solved using secant iteration, which works with any saturation vapor pressure formulation.

If the air is saturated or supersaturated ($pᵛ ≥ pᵛ⁺(T)$), the dewpoint equals the actual temperature and T is returned.

Arguments

• pᵛ: Vapor pressure (Pa)
• T: Actual temperature (K), used as upper bound and first guess
• constants: ThermodynamicConstants
• surface: Surface type for saturation vapor pressure calculation

Keyword arguments

• tolerance: Relative tolerance for convergence (default: 1e-4)
• maxiter: Maximum number of iterations (default: 10)

dewpoint_temperature(
    pᵛ,
    T,
    constants,
    equilibrium::Breeze.Thermodynamics.AbstractPhaseEquilibrium;
    tolerance,
    maxiter
) -> Any

Compute the dewpoint temperature using a phase equilibrium model to determine the condensation surface based on temperature T.

equilibrated_surface(phase_equilibrium::AbstractPhaseEquilibrium, T)

Return the appropriate surface type for computing saturation vapor pressure given the phase equilibrium model and temperature T.

equilibrium_saturation_specific_humidity(
    T,
    p,
    qᵗ,
    constants,
    surface
) -> Any

Compute the equilibrium saturation specific humidity $qᵛ⁺$ for air at temperature T, reference pressure pᵣ, and total specific moisture qᵗ, over a given surface.

This function returns the correct saturation specific humidity in both saturated and unsaturated conditions:

• In unsaturated conditions ($qᵗ < qᵛ⁺$), all moisture is vapor and the density is computed assuming $qᵛ = qᵗ$.

• In saturated conditions ($qᵗ ≥ qᵛ⁺$), the vapor specific humidity equals the saturation value and the density is computed assuming $qᵛ = qᵛ⁺$.

The saturated formula corresponds to equation (37) in Pressel et al. (2015).

equilibrium_saturation_specific_humidity(
    T,
    pᵣ,
    qᵗ,
    constants,
    equilibrium::Breeze.Thermodynamics.AbstractPhaseEquilibrium
) -> Any

Compute the equilibrium saturation specific humidity using a phase equilibrium model to determine the condensation surface based on temperature T.

ice_latent_heat(T, constants::ThermodynamicConstants) -> Any

Return the latent heat of sublimation (vapor → ice) at temperature T.

The latent heat varies linearly with temperature:

\[ℒⁱ(T) = ℒⁱᵣ + (cᵖᵛ - cⁱ)(T - Tᵣ)\]

where $ℒⁱᵣ$ is the reference latent heat at the energy reference temperature $Tᵣ$, $cᵖᵛ$ is the heat capacity of vapor, and $cⁱ$ is the heat capacity of ice.

liquid_latent_heat(
    T,
    constants::ThermodynamicConstants
) -> Any

Return the latent heat of vaporization (vapor → liquid) at temperature T.

The latent heat varies linearly with temperature:

\[ℒˡ(T) = ℒˡᵣ + (cᵖᵛ - cˡ)(T - Tᵣ)\]

where $ℒˡᵣ$ is the reference latent heat at the energy reference temperature $Tᵣ$, $cᵖᵛ$ is the heat capacity of vapor, and $cˡ$ is the heat capacity of liquid water.

mixture_gas_constant(
    q::Breeze.Thermodynamics.MoistureMassFractions,
    constants::ThermodynamicConstants
) -> Any

Return the gas constant of moist air mixture [in J/(kg K)] given the specific humidity q and thermodynamic parameters constants.

The mixture gas constant is calculated as a weighted average of the dry air and water vapor gas constants:

\[Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ ,\]

where:

• Rᵈ is the dry air gas constant,
• Rᵛ is the water vapor gas constant,
• qᵈ is the mass fraction of dry air, and
• qᵛ is the mass fraction of water vapor.

Arguments

• q: the moisture mass fractions (vapor, liquid, and ice)
• constants: ThermodynamicConstants instance containing gas constants

mixture_gas_constant(
    r::Breeze.Thermodynamics.MoistureMixingRatio,
    constants::ThermodynamicConstants
) -> Any

Compute the gas constant of a moist air mixture given moisture mixing ratios.

Converts mixing ratios to mass fractions and calls mixture_gas_constant(q::MMF, constants).

mixture_heat_capacity(
    q::Breeze.Thermodynamics.MoistureMassFractions,
    constants::ThermodynamicConstants
) -> Any

Compute the heat capacity of a mixture of dry air, vapor, liquid, and ice, where the mass fractions of vapor, liquid, and ice are given by q. The heat capacity of moist air is the weighted sum of its constituents:

\[cᵖᵐ = qᵈ cᵖᵈ + qᵛ cᵖᵛ + qˡ cˡ + qⁱ cⁱ ,\]

where qᵛ = q.vapor, qˡ = q.liquid, qⁱ = q.ice are the mass fractions of vapor, liquid, and ice constituents, respectively, and qᵈ = 1 - qᵛ - qˡ - qⁱ is the mass fraction of dry air. The heat capacities cᵖᵈ, cᵖᵛ, cˡ, cⁱ are the heat capacities of dry air, vapor, liquid, and ice at constant pressure, respectively. The liquid and ice phases are assumed to be incompressible.

mixture_heat_capacity(
    r::Breeze.Thermodynamics.MoistureMixingRatio,
    constants::ThermodynamicConstants
) -> Any

Compute the heat capacity of a moist air mixture given moisture mixing ratios.

Converts mixing ratios to mass fractions and calls mixture_heat_capacity(q::MMF, constants).

relative_humidity(T, ρ, qᵛ, constants) -> Any
relative_humidity(T, ρ, qᵛ, constants, surface) -> Any

Compute the relative humidity as the ratio of vapor pressure to saturation vapor pressure:

\[ℋ = pᵛ / pᵛ⁺ = qᵛ / qᵛ⁺\]

saturation_specific_humidity(
    T,
    ρ,
    constants,
    surface
) -> Any

Compute the saturation specific humidity for a gas at temperature T, total density ρ, constantsdynamics, and over surface via:

\[qᵛ⁺ = pᵛ⁺ / (ρ Rᵛ T) ,\]

where $pᵛ⁺$ is the saturation_vapor_pressure over surface, $ρ$ is total density, and $Rᵛ$ is the specific gas constant for water vapor.

Examples

First we compute the saturation specific humidity over a liquid surface:

using Breeze
using Breeze.Thermodynamics: PlanarLiquidSurface, PlanarIceSurface, PlanarMixedPhaseSurface

constants = ThermodynamicConstants()
T = 288.0 # Room temperature (K)
p = 101325.0 # Mean sea-level pressure
Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(constants)
q = zero(Breeze.Thermodynamics.MoistureMassFractions{Float64})
ρ = Breeze.Thermodynamics.density(T, p, q, constants)
qᵛ⁺ˡ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, constants, PlanarLiquidSurface())

# output
0.010359995391195264

Note, this is slightly smaller than the saturation specific humidity over an ice surface:

julia> qᵛ⁺ˡ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, constants, PlanarIceSurface())
0.011945100768555072

If a medium contains a mixture of 40% water and 60% ice that has (somehow) acquired thermodynamic equilibrium, we can compute the saturation specific humidity over the mixed phase surface,

mixed_surface = PlanarMixedPhaseSurface(0.4)
qᵛ⁺ᵐ = Breeze.Thermodynamics.saturation_specific_humidity(T, ρ, constants, mixed_surface)

# output
0.01128386068542303

saturation_vapor_pressure(
    T,
    constants::Breeze.Thermodynamics.ClausiusClapeyronThermodynamicConstants,
    surface
) -> Any

Compute the saturation vapor pressure $pᵛ⁺$ over a surface labeled $β$ (for example, a planar liquid surface, or curved ice surface) using the Clausius-Clapeyron relation,

\[𝖽pᵛ⁺ / 𝖽T = pᵛ⁺ ℒᵝ(T) / (Rᵛ T^2) ,\]

where the temperature-dependent latent heat of the surface is $ℒᵝ(T)$.

Using a model for the latent heat that is linear in temperature, eg

\[ℒᵝ = ℒᵝ₀ + Δcᵝ T,\]

where $ℒᵝ₀ ≡ ℒᵝ(T=0)$ is the latent heat at absolute zero and $Δcᵝ ≡ cᵖᵛ - cᵝ$ is the constant difference between the vapor specific heat and the specific heat of phase $β$.

Note that we typically parameterize the latent heat in terms of a reference temperature $T = Tᵣ$ that is well above absolute zero. In that case, the latent heat is written

\[ℒᵝ = ℒᵝᵣ + Δcᵝ (T - Tᵣ) \qquad \text{and} \qquad ℒᵝ₀ = ℒᵝᵣ - Δcᵝ Tᵣ .\]

Integrating the Clausius-Clapeyron relation with a temperature-linear latent heat model, from the triple point pressure and temperature $(pᵗʳ, Tᵗʳ)$ to pressure $pᵛ⁺$ and temperature $T$, we obtain

\[\log(pᵛ⁺ / pᵗʳ) = - ℒᵝ₀ / (Rᵛ T) + ℒᵝ₀ / (Rᵛ Tᵗʳ) + \log \left[ (Δcᵝ / Rᵛ) (T / Tᵗʳ) \right] ,\]

which then becomes

\[pᵛ⁺(T) = pᵗʳ (T / Tᵗʳ)^{Δcᵝ / Rᵛ} \exp \left [ (1/Tᵗʳ - 1/T) ℒᵝ₀ / Rᵛ \right ] .\]

saturation_vapor_pressure(
    T,
    constants::Breeze.Thermodynamics.TetensFormulaThermodynamicConstants,
    surface::Breeze.Thermodynamics.PlanarMixedPhaseSurface
) -> Any

Compute the saturation vapor pressure over a mixed-phase surface by linearly interpolating between liquid and ice saturation vapor pressures based on the liquid fraction.

saturation_vapor_pressure(
    T,
    constants::Breeze.Thermodynamics.TetensFormulaThermodynamicConstants,
    _::PlanarIceSurface
) -> Any

Compute the saturation vapor pressure over a planar ice surface using Tetens' empirical formula with ice coefficients from Murray (1967):

\[pᵛ⁺(T) = pᵛ⁺ᵣ \exp \left( aⁱ \frac{T - Tᵣ}{T - δTⁱ} \right)\]

References

• Tetens, O. (1930). Über einige meteorologische Begriffe. Zeitschrift für Geophysik 6, 297–309.
• Murray, F. W. (1967). On the computation of saturation vapor pressure. Journal of Applied Meteorology 6, 203–204.

saturation_vapor_pressure(
    T,
    constants::Breeze.Thermodynamics.TetensFormulaThermodynamicConstants,
    _::PlanarLiquidSurface
) -> Any

Compute the saturation vapor pressure over a planar liquid surface using Tetens' empirical formula:

\[pᵛ⁺(T) = pᵛ⁺ᵣ \exp \left( aˡ \frac{T - Tᵣ}{T - δTˡ} \right)\]

supersaturation(
    T,
    ρ,
    q::Breeze.Thermodynamics.MoistureMassFractions,
    constants,
    surface
) -> Any

Compute the supersaturation $𝒮 = pᵛ/pᵛ⁺ - 1$ over a given surface.

• $𝒮 < 0$ indicates subsaturation (evaporation conditions)
• $𝒮 = 0$ indicates saturation (equilibrium)
• $𝒮 > 0$ indicates supersaturation (condensation conditions)

Arguments

• T: Temperature
• ρ: Total air density
• q: MoistureMassFractions containing vapor, liquid, and ice mass fractions
• constants: ThermodynamicConstants
• surface: Surface type (e.g., PlanarLiquidSurface(), PlanarIceSurface())

surface_density(p₀, θ₀, pˢᵗ, constants)

Compute the surface air density from surface pressure p₀, potential temperature θ₀, standard pressure pˢᵗ, and thermodynamic constants using the ideal gas law for dry air.

The temperature is computed from potential temperature using the Exner function: T₀ = Π₀ * θ₀ where Π₀ = (p₀ / pˢᵗ)^(Rᵈ/cᵖᵈ).

surface_density(p₀, T₀, constants)

Compute the surface air density from surface pressure p₀, surface temperature T₀, and thermodynamic constants using the ideal gas law for dry air.

surface_density(reference_state)

Return the density at z=0 by interpolating the reference density field to the surface.

temperature_from_potential_temperature(θ, p, constants; pˢᵗ=1e5, qᵛ=0)

Compute temperature from potential temperature and pressure.

This is a convenience function that constructs a LiquidIcePotentialTemperatureState with no condensate and computes temperature using the standard thermodynamic relations.

Arguments

• θ: Potential temperature [K]
• p: Pressure [Pa]
• constants: Thermodynamic constants

Keyword Arguments

• pˢᵗ: Standard pressure for potential temperature definition [Pa] (default: 1e5)
• qᵛ: Specific humidity [kg/kg] (default: 0, dry air)

vapor_pressure(T, ρ, qᵛ, constants) -> Any

Compute the vapor pressure from the ideal gas law:

\[pᵛ = ρ qᵛ Rᵛ T\]

TimeSteppers module for Breeze.jl

Provides time stepping schemes for AtmosphereModel, including:

• SSPRungeKutta3: Standard SSP RK3 scheme for explicit time stepping
• AcousticSSPRungeKutta3: SSP RK3 with acoustic substepping for compressible dynamics
• AcousticRungeKutta3: Wicker-Skamarock RK3 with acoustic substepping for compressible dynamics

struct AcousticRungeKutta3{FT, U0, TG, TI, AS} <: Oceananigans.TimeSteppers.AbstractTimeStepper

Wicker-Skamarock third-order Runge-Kutta time stepper with acoustic substepping for fully compressible dynamics.

Unlike AcousticSSPRungeKutta3 which uses convex combinations, this scheme uses stage fractions Δt/3, Δt/2, Δt:

• Stage 1: $U^* = U^n + (Δt/3) \, R(U^n)$
• Stage 2: $U^{**} = U^n + (Δt/2) \, R(U^*)$
• Stage 3: $U^{n+1} = U^n + Δt \, R(U^{**})$

Each stage evaluates the RHS at the current stage state, then resets to $U^n$ and advances by $β Δt$. The absence of convex combinations makes this scheme compatible with split-explicit acoustic substepping, allowing the full pressure gradient and buoyancy to be included in the slow tendency.

This is the scheme used by WRF and CM1 for compressible atmospheric dynamics.

Fields

• β₁, β₂, β₃: Stage fractions (1/3, 1/2, 1)
• U⁰: Storage for state at beginning of time step
• Gⁿ: Tendency fields at current stage
• implicit_solver: Optional implicit solver for diffusion
• substepper: AcousticSubstepper for acoustic substepping infrastructure

References

Wicker, L.J. and Skamarock, W.C. (2002). Time-Splitting Methods for Elastic Models Using Forward Time Schemes. Monthly Weather Review, 130, 2088-2097.

Klemp, J.B., Skamarock, W.C. and Dudhia, J. (2007). Conservative Split-Explicit Time Integration Methods for the Compressible Nonhydrostatic Equations. Monthly Weather Review, 135, 2897-2913.

AcousticRungeKutta3(grid, prognostic_fields;
                     dynamics,
                     implicit_solver = nothing,
                     Gⁿ = map(similar, prognostic_fields))

Construct an AcousticRungeKutta3 time stepper for fully compressible dynamics.

Keyword Arguments

• dynamics: The CompressibleDynamics object containing the time_discretization.
• implicit_solver: Optional implicit solver for diffusion. Default: nothing
• Gⁿ: Tendency fields at current stage. Default: similar to prognostic_fields

struct AcousticSSPRungeKutta3{FT, U0, TG, TI, AS} <: Oceananigans.TimeSteppers.AbstractTimeStepper

A strong stability preserving (SSP) third-order Runge-Kutta time stepper with acoustic substepping for fully compressible dynamics.

This time stepper implements the Wicker-Skamarock scheme used in CM1:

• Outer RK3 loop for slow tendencies (advection, buoyancy, turbulence)
• Inner acoustic substep loop for fast tendencies (pressure gradient, compression)

The acoustic substepping separates time scales:

• Slow modes (advection, buoyancy): CFL ≈ 10-20 m/s → Δtˢˡᵒʷ ~ 1-10 s
• Fast modes (acoustic): CFL ≈ 340 m/s → Δtˢ ~ 0.1-0.3 s

By substepping the fast modes, we can use ~6 acoustic substeps per slow step instead of reducing the overall time step by a factor of ~30.

Fields

• α¹, α², α³: SSP RK3 stage coefficients (1, 1/4, 2/3)
• U⁰: Storage for state at beginning of time step
• Gⁿ: Tendency fields at current stage
• implicit_solver: Optional implicit solver for diffusion
• substepper: AcousticSubstepper for acoustic substepping infrastructure

AcousticSSPRungeKutta3(grid, prognostic_fields;
                       dynamics,
                       implicit_solver = nothing,
                       Gⁿ = map(similar, prognostic_fields))

Construct an AcousticSSPRungeKutta3 time stepper for fully compressible dynamics.

This combines the SSP RK3 scheme from Shu and Osher (1988) with acoustic substepping from Wicker and Skamarock (2002).

The acoustic substepping parameters are configured via the time_discretization field of the CompressibleDynamics object passed as dynamics.

Keyword Arguments

• dynamics: The CompressibleDynamics object containing the time_discretization.
• implicit_solver: Optional implicit solver for diffusion. Default: nothing
• Gⁿ: Tendency fields at current stage. Default: similar to prognostic_fields

References

Shu, C.-W., & Osher, S. (1988). Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77(2), 439-471.

Wicker, L.J. and Skamarock, W.C. (2002). Time-Splitting Methods for Elastic Models Using Forward Time Schemes. Monthly Weather Review, 130, 2088-2097.

struct SSPRungeKutta3{FT, U0, TG, TI} <: Oceananigans.TimeSteppers.AbstractTimeStepper

A strong stability preserving (SSP) third-order Runge-Kutta time stepper.

This time stepper uses the classic SSP RK3 scheme (Shu-Osher 2006 form):

\[\begin{align*} u^{(1)} &= u^{(0)} + Δt L(u^{(0)}) \\ u^{(2)} &= \frac{3}{4} u^{(0)} + \frac{1}{4} u^{(1)} + \frac{1}{4} Δt L(u^{(1)}) \\ u^{(3)} &= \frac{1}{3} u^{(0)} + \frac{2}{3} u^{(2)} + \frac{2}{3} Δt L(u^{(2)}) \end{align*}\]

where $L$ above is the right-hand-side, e.g., $\partial_t u = L(u)$.

Each stage can be written in the form:

\[u^{(m)} = (1 - α) u^{(0)} + α [u^{(m-1)} + Δt L(u^{(m-1)})]\]

with $α = 1, 1/4, 2/3$ for stages 1, 2, 3 respectively.

This scheme has CFL coefficient = 1 and is TVD (total variation diminishing).

Fields

• α¹, α², α³: Stage coefficients (1, 1/4, 2/3)
• U⁰: Storage for state at beginning of time step
• Gⁿ: Tendency fields at current stage
• implicit_solver: Optional implicit solver for diffusion

SSPRungeKutta3(grid, prognostic_fields;
               implicit_solver = nothing,
               Gⁿ = map(similar, prognostic_fields))

Construct an SSPRungeKutta3 on grid with prognostic_fields as described by Shu and Osher (1988).

Keyword Arguments

• implicit_solver: Optional implicit solver for diffusion. Default: nothing
• Gⁿ: Tendency fields at current stage. Default: similar to prognostic_fields

References

Shu, C.-W., & Osher, S. (1988). Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77(2), 439-471.

maybe_initialize_state!(model, callbacks)

Update the model state at iteration 0, in case run! is not used.

ssp_rk3_substep!(model, Δt, α)

Apply an SSP RK3 substep with coefficient α:

u^(m) = (1 - α) * u^(0) + α * (u^(m-1) + Δt * G)

where u^(0) is stored in the time stepper, u^(m-1) is the current field value, and G is the current tendency.

store_initial_state!(model)

Copy prognostic fields to U⁰ storage for use in later RK3 stages.

PiecewiseStretchedDiscretization(; z, Δz)

Construct a stretched vertical grid where the spacing varies piecewise-linearly between breakpoints. The grid spacing is specified at breakpoint heights z, and linearly interpolated between them.

Between breakpoints where Δz values are equal, the grid is uniform. Where they differ, the spacing transitions linearly.

The result behaves as a vector of face positions and can be passed directly to RectilinearGrid as a coordinate argument.

Keyword Arguments

• z: sorted vector of breakpoint heights (length ≥ 2)
• Δz: vector of grid spacings at each breakpoint (same length as z, all positive)

Examples

A three-region grid with uniform fine spacing, a linear transition, and uniform coarse spacing (as used for tropical cyclone simulations):

z = PiecewiseStretchedDiscretization(
    z  = [0, 1000, 3500, 28000],
    Δz = [62.5, 62.5, 2000, 2000]
)
Nz = length(z) - 1
grid = RectilinearGrid(arch; size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z)

A four-region grid for deep convection (fine near surface, transition to moderate, uniform through the troposphere, then stretched to the model top):

z = PiecewiseStretchedDiscretization(
    z  = [0, 1275, 5100, 18000, 27000],
    Δz = [50, 50, 100, 100, 300]
)

buoyancy_forceᶜᶜᶜ(
    i,
    j,
    k,
    grid,
    dynamics::AnelasticDynamics,
    temperature,
    specific_moisture,
    microphysics,
    microphysical_fields,
    constants
) -> Any

Compute the buoyancy force density for anelastic dynamics at cell center (i, j, k).

The anelastic buoyancy force is the gravitational force on the density anomaly:

\[-g \rho' = -g (\rho - \rho_r)\]

where $\rho = p_r / (R^m T)$ is the in-situ density from the ideal gas law, and $\rho_r = p_r / (R^m_r T_r)$ is the reference density. Substituting:

\[\rho' = \frac{p_r}{R^m T} - \frac{p_r}{R^m_r T_r} = \frac{p_r}{R^m_r T_r} \left( \frac{R^m_r T_r}{R^m T} - 1 \right) = \rho_r \left( \frac{R^m_r T_r}{R^m T} - 1 \right)\]

This "perturbation form" avoids subtracting two large, nearly-equal numbers ($p_r / (R^m T) - \rho_r$), which causes catastrophic cancellation when $T \approx T_r$. Instead, the ratio $R^m_r T_r / (R^m T)$ is close to 1, and the subtraction of 1 preserves relative precision.

Here $R^m = q^d R^d + q^v R^v$ is the mixture gas constant for the current moisture state and $R^m_r$ is the mixture gas constant for the reference moisture state.

compute_pressure_correction!(
    model::AtmosphereModel{<:AnelasticDynamics},
    Δt
)

Compute the pressure correction for anelastic dynamics by solving the pressure Poisson equation.

default_dynamics(
    grid,
    constants
) -> AnelasticDynamics{R, Nothing} where R<:(ReferenceState{_A, P, D, T, QV, QL, QI} where {_A, P<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), D<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), T<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), QV<:(Oceananigans.Fields.ZeroField{T, 3} where T), QL<:(Oceananigans.Fields.ZeroField{T, 3} where T), QI<:(Oceananigans.Fields.ZeroField{T, 3} where T)})

Construct a "stub" AnelasticDynamics with just the reference_state. The pressure anomaly field is materialized later in the model constructor.

dynamics_density(dynamics::AnelasticDynamics) -> Any

Return the reference density field for AnelasticDynamics.

For anelastic models, the dynamics density is the time-independent reference state density ρᵣ(z).

dynamics_pressure(dynamics::AnelasticDynamics) -> Any

Return the reference pressure field for AnelasticDynamics.

For anelastic models, the dynamics pressure is the time-independent hydrostatic reference state pressure pᵣ(z).

initialize_model_thermodynamics!(
    model::AtmosphereModel{<:AnelasticDynamics}
)

Initialize thermodynamic state for anelastic models. Sets the initial potential temperature to the reference state value.

make_pressure_correction!(
    model::AtmosphereModel{<:AnelasticDynamics},
    Δt
)

Update the predictor momentum $(ρu, ρv, ρw)$ with the non-hydrostatic pressure via

\[(\rho\boldsymbol{u})^{n+1} = (\rho\boldsymbol{u})^n - \Delta t \, \rho_r \boldsymbol{\nabla} \left( \alpha_r p_{nh} \right)\]

materialize_dynamics(
    dynamics::AnelasticDynamics,
    grid,
    boundary_conditions,
    thermodynamic_constants
) -> AnelasticDynamics{_A, P} where {_A, P<:(Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B})}

Materialize a stub AnelasticDynamics into a full dynamics object with the pressure anomaly field.

mean_pressure(dynamics::AnelasticDynamics) -> Any

Return the mean (reference) pressure field for AnelasticDynamics, in Pa.

pressure_anomaly(dynamics::AnelasticDynamics) -> Any

Return the non-hydrostatic pressure anomaly for AnelasticDynamics, in Pa. Note: the internal field stores the kinematic pressure p'/ρᵣ; this function returns ρᵣ * p'/ρᵣ = p' in Pa.

standard_pressure(dynamics::AnelasticDynamics) -> Any

Return the standard pressure from the reference state for potential temperature calculations.

surface_pressure(dynamics::AnelasticDynamics) -> Any

Return the surface pressure from the reference state for boundary condition regularization.

total_pressure(dynamics::AnelasticDynamics) -> Any

Return the total pressure for AnelasticDynamics, in Pa. This is p = p̄ + p', where p̄ is the hydrostatic reference pressure and p' is the non-hydrostatic pressure anomaly.

abstract type AbstractOptics

Abstract type representing optics for RadiativeTransferModel.

additional_dynamics_field_names(dynamics)

Return a tuple of additional (diagnostic) field names for the dynamics.

additional_thermodynamic_field_names(formulation)

Return a tuple of additional (diagnostic) field names for the given thermodynamic formulation. Accepts a Symbol, Val(Symbol), or formulation struct.

adjust_thermodynamic_state(
    state,
    scheme::Nothing,
    thermo
) -> Any

Adjust the thermodynamic state according to the scheme. For example, if scheme isa SaturationAdjustment, then this function will adjust and return a new thermodynamic state given the specifications of the saturation adjustment scheme.

If a scheme is non-adjusting, we just return state.

buoyancy_forceᶜᶜᶠ(i, j, k, grid, args...) -> Any

Interpolate buoyancy force to z-face location.

cloud_ice_effective_radius(
    i,
    j,
    k,
    grid,
    effective_radius_model::ConstantRadiusParticles,
    args...
) -> Any

Return the effective radius of cloud ice particles in microns (μm).

This function dispatches on the effective_radius_model argument. The default implementation for ConstantRadiusParticles returns a constant value.

Microphysics schemes can extend this function to provide diagnosed effective radii based on cloud properties.

cloud_liquid_effective_radius(
    i,
    j,
    k,
    grid,
    effective_radius_model::ConstantRadiusParticles,
    args...
) -> Any

Return the effective radius of cloud liquid droplets in microns (μm).

This function dispatches on the effective_radius_model argument. The default implementation for ConstantRadiusParticles returns a constant value.

Microphysics schemes can extend this function to provide diagnosed effective radii based on cloud properties.

collect_prognostic_fields(formulation, dynamics, momentum, moisture_density, microphysical_fields, tracers)

Collect all prognostic fields into a single NamedTuple.

compute_auxiliary_dynamics_variables!(model)

Compute auxiliary (diagnostic) variables specific to the dynamics formulation.

For anelastic dynamics, this is a no-op (pressure is computed during time-stepping via the pressure Poisson equation).

For compressible dynamics, this computes the pressure field from the equation of state:

\[p = ρ R^m T\]

where $R^m$ is the mixture gas constant.

compute_auxiliary_thermodynamic_variables!(formulation, dynamics, i, j, k, grid)

Compute auxiliary thermodynamic variables from prognostic fields at grid point (i, j, k).

compute_auxiliary_variables!(model)

Compute auxiliary model variables:

• velocities from momentum and density (eg $u = ρu / ρ$)

• thermodynamic variables from the prognostic thermodynamic state,

  • temperature $T$, possibly involving saturation adjustment
  • specific thermodynamic variable ($e = ρe / ρ$ or $θ = ρθ / ρ$)
  • moisture mass fraction $qᵗ = ρqᵗ / ρ$

compute_dynamics_tendency!(model)

Compute tendencies for dynamics-specific prognostic fields.

For anelastic dynamics, this is a no-op (no prognostic density). For compressible dynamics, this computes the density tendency from the continuity equation:

\[\partial_t \rho = -\boldsymbol{\nabla \cdot} (\rho \boldsymbol{u})\]

compute_forcings!(model)

Compute forcing-specific quantities needed before tendency calculation. For example, SubsidenceForcing requires horizontal averages of the fields being advected.

compute_thermodynamic_tendency!(model, common_args)

Compute the thermodynamic tendency. Dispatches on the thermodynamic formulation type.

compute_velocities!(model::AtmosphereModel)

Compute velocities from momentum: u = ρu / ρ for each velocity component.

default_dynamics(grid, constants)

Return the default dynamics for the given grid and thermodynamic constants.

default_timestepper(dynamics) -> Symbol

Return the default timestepper symbol for the given dynamics.

For anelastic dynamics, returns :SSPRungeKutta3. For compressible dynamics, returns :AcousticSSPRungeKutta3 (acoustic substepping).

diagnose_thermodynamic_state(i, j, k, grid, formulation, dynamics, q)

Diagnose the thermodynamic state at grid point (i, j, k) from the given formulation, dynamics, and pre-computed moisture mass fractions q.

Note: This function does NOT compute moisture fractions internally to avoid circular dependencies. The caller is responsible for computing q = grid_moisture_fractions(...) before passing q to this function.

dynamics_pressure_solver(dynamics, grid)

Create the pressure solver for the given dynamics. Returns nothing for dynamics that do not require a pressure solver (e.g., compressible).

dynamics_prognostic_fields(dynamics)

Return a NamedTuple of prognostic fields specific to the dynamics formulation.

For anelastic dynamics, returns an empty NamedTuple. For compressible dynamics, returns (ρ=density_field,).

extract_microphysical_prognostics(
    i,
    j,
    k,
    microphysics,
    μ_fields
) -> NamedTuple

Extract prognostic microphysical variables at grid point (i, j, k) into a NamedTuple of scalar values.

Uses prognostic_field_names to determine which fields to extract. The result is a NamedTuple with density-weighted values (e.g., (ρqᶜˡ=..., ρqʳ=...)).

This function enables a generic grid-indexed microphysical_state that extracts prognostics and delegates to the gridless version.

grid_microphysical_state(i, j, k, grid, microphysics, μ_fields, ρ, 𝒰, velocities)

Build an AbstractMicrophysicalState (ℳ) at grid point (i, j, k).

This is the grid-indexed wrapper that:

1. Extracts prognostic values from μ_fields via extract_microphysical_prognostics
2. Calls the gridless microphysical_state(microphysics, ρ, μ, 𝒰, velocities)

Microphysics schemes should implement the gridless version, not this one.

Arguments

• i, j, k: Grid indices
• grid: The computational grid
• microphysics: The microphysics scheme
• μ_fields: NamedTuple of microphysical fields
• ρ: Local density (scalar)
• 𝒰: Thermodynamic state
• velocities: Velocity fields (u, v, w). Velocities are interpolated to cell centers for use by microphysics schemes (e.g., aerosol activation uses vertical velocity).

Returns

An AbstractMicrophysicalState subtype containing the local microphysical variables.

See also microphysical_tendency, AbstractMicrophysicalState.

grid_moisture_fractions(
    i,
    j,
    k,
    grid,
    microphysics,
    ρ,
    qᵗ,
    μ_fields
) -> Breeze.Thermodynamics.MoistureMassFractions

Grid-indexed version of moisture_fractions.

This is the generic wrapper that:

1. Extracts prognostic values from μ_fields via extract_microphysical_prognostics
2. Builds the microphysical state via microphysical_state with 𝒰 = nothing
3. Calls moisture_fractions

This works for non-equilibrium schemes where cloud condensate is prognostic. Non-equilibrium schemes don't need 𝒰 to build their state (they use prognostic fields).

Saturation adjustment schemes should override this to read from diagnostic fields.

initialize_model_microphysical_fields!(fields, _::Nothing)

Initialize default values for microphysical fields after materialization.

Sets ρnᵃ (aerosol number density) to initial_aerosol_number(microphysics) if the field exists. All other microphysical fields remain at zero. Users can override with set!.

initialize_model_thermodynamics!(model)

Initialize the thermodynamic state for a newly constructed model. For anelastic dynamics, sets initial θ to the reference potential temperature. For compressible dynamics, no default initialization is performed.

materialize_dynamics(dynamics_stub, grid, boundary_conditions, thermodynamic_constants, microphysics=nothing)

Materialize a dynamics stub into a complete dynamics object with all required fields.

The microphysics argument is optional and used by dynamics types that need to know the microphysics scheme to create appropriate prognostic state (e.g., ParcelDynamics).

materialize_formulation(formulation, dynamics, grid, boundary_conditions)

Materialize a thermodynamic formulation from a Symbol (or formulation struct) into a complete formulation with all required fields.

Valid symbols:

• :LiquidIcePotentialTemperature, :θ, :ρθ, :PotentialTemperature → LiquidIcePotentialTemperatureFormulation
• :StaticEnergy, :e, :ρe → StaticEnergyFormulation

materialize_microphysical_fields(
    microphysics::Nothing,
    grid,
    boundary_conditions
) -> @NamedTuple{}

Build microphysical fields associated with microphysics on grid and with user defined boundary_conditions.

materialize_momentum_and_velocities(dynamics, grid, boundary_conditions)

Create momentum and velocity fields for the given dynamics.

materialize_velocities(velocities, grid)

Create velocity fields from a velocity specification (e.g., PrescribedVelocityFields).

maybe_adjust_thermodynamic_state(
    state,
    _::Nothing,
    qᵗ,
    constants
) -> Any

Possibly apply saturation adjustment. If a microphysics scheme does not invoke saturation adjustment, just return the state unmodified.

This function takes the thermodynamic state, microphysics scheme, total moisture, and thermodynamic constants. Schemes that use saturation adjustment override this to adjust the moisture partition. Non-equilibrium schemes simply return the state unchanged.

microphysical_velocities(
    microphysics::Nothing,
    microphysical_fields,
    name
)

Return the microphysical velocities associated with microphysics, microphysical_fields, and tracer name.

Must be either nothing, or a NamedTuple with three components u, v, w. The velocities are added to the bulk flow velocities for advecting the tracer. For example, the terminal velocity of falling rain.

microphysics_model_update!(microphysics::Nothing, model)

Apply microphysics model update for the given microphysics scheme.

This function is called during update_state! to apply microphysics processes that operate on the full model state (not the tendency fields). Specific microphysics schemes should extend this function.

moisture_fractions(
    _::Nothing,
    ℳ,
    qᵗ
) -> Breeze.Thermodynamics.MoistureMassFractions

Compute MoistureMassFractions from a microphysical state ℳ and total moisture qᵗ.

This is the state-based (gridless) interface for computing moisture fractions. Microphysics schemes should extend this method to partition moisture based on their prognostic variables.

The default implementation for Nothing microphysics assumes all moisture is vapor.

prognostic_dynamics_field_names(dynamics)

Return a tuple of prognostic field names specific to the dynamics formulation.

For anelastic dynamics, returns an empty tuple (no prognostic density). For compressible dynamics, returns (:ρ,) for prognostic density.

prognostic_field_names(_::Nothing) -> Tuple{}

Return tuple() - zero-moment scheme has no prognostic variables.

prognostic_momentum_field_names(dynamics)

Return a tuple of prognostic momentum field names.

For prognostic dynamics (anelastic, compressible), returns (:ρu, :ρv, :ρw). For kinematic dynamics (prescribed velocities), returns an empty tuple.

prognostic_thermodynamic_field_names(formulation)

Return a tuple of prognostic field names for the given thermodynamic formulation. Accepts a Symbol, Val(Symbol), or formulation struct.

set_momentum!(model, name, value)

Set the momentum component name (:ρu, :ρv, or :ρw) to value.

set_thermodynamic_variable!(model, variable_name, value)

Set a thermodynamic variable (e.g., :θ, :T, :e, :ρθ, :ρe) from the given value. Dispatches on the thermodynamic formulation type and variable name.

set_velocity!(model, name, value)

Set the velocity component name (:u, :v, or :w) to value. Also updates the corresponding momentum field.

settable_specific_microphysical_names(
    microphysics
) -> Tuple{}

Return a tuple of specific (non-density-weighted) names that can be set for the given microphysics scheme. These are derived from the prognostic field names by removing the 'ρ' prefix.

For mass fields (e.g., ρqᶜˡ → qᶜˡ) and number fields (e.g., ρnᶜˡ → nᶜˡ).

specific_to_density_weighted(
    name::Symbol
) -> Union{Nothing, Symbol}

Convert a specific microphysical variable name to its density-weighted counterpart. For example, :qᶜˡ → :ρqᶜˡ, :qʳ → :ρqʳ, :nᶜˡ → :ρnᶜˡ.

Returns nothing if the name doesn't start with 'q' or 'n'.

standard_pressure(dynamics)

Return the standard pressure used for potential temperature calculations. Default is 100000 Pa (1000 hPa).

surface_pressure(dynamics)

Return the surface pressure used for boundary condition regularization. For anelastic dynamics, this is the reference state surface pressure. For compressible dynamics, this may be a constant or computed value.

update_dynamics_with_velocities(dynamics, velocities)

Update dynamics with velocity specification. Default is a no-op. For PrescribedDynamics, stores the PrescribedVelocityFields for dispatch.

update_radiation!(radiation, model)

Update the radiative fluxes from the current model state.

This is a stub function that does nothing by default. It is extended by radiation extensions (e.g., BreezeRRTMGPExt) to compute radiative transfer.

validate_velocity_boundary_conditions(dynamics, user_boundary_conditions)

Validate that velocity boundary conditions are only provided for dynamics that support them.

By default, throws an error if the user provides boundary conditions for :u, :v, or :w, since velocity is a diagnostic field for most dynamics (e.g., anelastic, compressible).

For PrescribedDynamics, velocity boundary conditions are allowed since velocities are regular fields that can be set directly.

velocity_boundary_condition_names(dynamics)

Return a tuple of velocity field names that need default boundary conditions.

For most dynamics (anelastic, compressible), velocities are diagnostic and their boundary conditions are created internally. Returns an empty tuple.

For PrescribedDynamics, velocities are regular fields that can have user-provided boundary conditions, so this returns (:u, :v, :w).

w_buoyancy_forceᶜᶜᶠ(i, j, k, grid, w, args...) -> Any

Compute the product of vertical velocity and buoyancy force at z-face location. Used for the buoyancy flux term in the energy equation.

x_pressure_gradient(i, j, k, grid, dynamics)

Return the x-component of the pressure gradient force at (Face, Center, Center).

For anelastic dynamics, returns zero (pressure is handled via projection). For compressible dynamics, returns -∂p/∂x.

y_pressure_gradient(i, j, k, grid, dynamics)

Return the y-component of the pressure gradient force at (Center, Face, Center).

For anelastic dynamics, returns zero (pressure is handled via projection). For compressible dynamics, returns -∂p/∂y.

z_pressure_gradient(i, j, k, grid, dynamics)

Return the z-component of the pressure gradient force at (Center, Center, Face).

For anelastic dynamics, returns zero (pressure is handled via projection). For compressible dynamics, returns -∂p/∂z.

∇_dot_Jᶜ(i, j, k, grid, ρ, closure::AbstractTurbulenceClosure, closure_fields,
         id, c, clock, model_fields, buoyancy)

Return the discrete divergence of the dynamic scalar flux Jᶜ = ρ jᶜ, where jᶜ is the "kinematic scalar flux", using area-weighted differences divided by cell volume. Similar to Oceananigans' ∇_dot_qᶜ signature with the additional density factor ρ, where in Oceananigans qᶜ is the kinematic tracer flux.

set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...)

Set variables in an AtmosphereModel.

Keyword Arguments

Variables are set via keyword arguments. Supported variables include:

Prognostic variables (density-weighted):

• ρu, ρv, ρw: momentum components
• ρqᵗ: total moisture density
• Prognostic microphysical variables
• Prognostic user-specified tracer fields

Settable thermodynamic variables:

• T: in-situ temperature
• θ: potential temperature
• θˡⁱ: liquid-ice potential temperature
• e: static energy
• ρθ: potential temperature density
• ρθˡⁱ: liquid-ice potential temperature density
• ρe: static energy density (for StaticEnergyThermodynamics)

Diagnostic variables (specific, i.e., per unit mass):

• u, v, w: velocity components (sets both velocity and momentum)
• qᵗ: total specific moisture (sets both specific and density-weighted moisture)
• ℋ: relative humidity (sets total moisture via qᵗ = ℋ * qᵛ⁺, where qᵛ⁺ is the saturation specific humidity at the current temperature). Relative humidity is in the range [0, 1]. For models with saturation adjustment microphysics, ℋ > 1 throws an error since the saturation adjustment would immediately reduce it to 1.

Specific microphysical variables (automatically converted to density-weighted):

• qᶜˡ: specific cloud liquid (sets ρqᶜˡ = ρᵣ * qᶜˡ)
• qʳ: specific rain (sets ρqʳ = ρᵣ * qʳ)
• nᶜˡ: specific cloud liquid number [1/kg] (sets ρnᶜˡ = ρᵣ * nᶜˡ)
• nʳ: specific rain number [1/kg] (sets ρnʳ = ρᵣ * nʳ)
• Other prognostic microphysical variables with the ρ prefix removed

Options

• enforce_mass_conservation: If true (default), applies a pressure correction to ensure the velocity field satisfies the anelastic continuity equation.

compute_flux_bc_tendencies!(model::AtmosphereModel)

Apply boundary conditions by adding flux divergences to the right-hand-side.

dynamics_density_for_potential_temperature(dynamics) -> Any

Return the density field used for computing potential temperature. For anelastic dynamics, this is the reference density.

dynamics_pressure_for_potential_temperature(dynamics) -> Any

Return the pressure field used for computing potential temperature. For anelastic dynamics, this is the reference pressure.

dynamics_standard_pressure(dynamics) -> Any

Return the standard pressure for potential temperature calculations.

materialize_atmosphere_model_boundary_conditions(
    boundary_conditions,
    grid,
    formulation,
    dynamics,
    microphysics,
    surface_pressure,
    thermodynamic_constants,
    microphysical_fields,
    specific_moisture,
    temperature
) -> NamedTuple

Regularize boundary conditions for AtmosphereModel. This function walks through all boundary conditions and calls materialize_atmosphere_boundary_condition on each one, allowing specialized handling for bulk flux boundary conditions and other atmosphere-specific boundary condition types.

If formulation is :LiquidIcePotentialTemperature and ρe boundary conditions are provided, they are automatically converted to ρθ boundary conditions using EnergyFluxBoundaryCondition.

bulk_richardson_number(h, θᵥ, θᵥ₀, U, U_min) -> Any
bulk_richardson_number(h, θᵥ, θᵥ₀, U, U_min, g) -> Any

Compute bulk Richardson number: Riᴮ = (g/θ̄ᵥ) × h × (θᵥ - θᵥ₀) / U²

Wind speed is clamped to U_min to avoid singularity.

Arguments

• h: Measurement height (m)
• θᵥ: Virtual potential temperature at measurement height (K)
• θᵥ₀: Virtual potential temperature at surface (K)
• U: Wind speed (m/s)
• U_min: Minimum wind speed (m/s)
• g: Gravitational acceleration (m/s², default: 9.81)

bulk_to_flux_richardson_number(Riᴮ, α, β, mapping) -> Any

Map bulk Richardson number $Riᴮ$ to the Monin-Obukhov stability parameter $ζ = z/L$ using the regression equations of Li et al. (2010).

Arguments

• Riᴮ: Bulk Richardson number
• α: $\ln(z / ℓ)$
• β: $\ln(ℓ / ℓ_h)$
• mapping: RichardsonNumberMapping with regression coefficients

integrated_stability_momentum(ζ, params) -> Any

Integrated stability function for momentum $Ψᴰ(ζ)$.

Note: $Ψᴰ$ corresponds to $Ψ_m$ in the literature.

integrated_stability_scalar(ζ, params) -> Any

Integrated stability function for scalars (heat, moisture) $Ψᵀ(ζ)$.

Note: $Ψᵀ$ corresponds to $Ψ_h$ in the literature.

neutral_coefficient_10m(polynomial, U₁₀, U_min) -> Any

Compute neutral transfer coefficient at 10 m using the Large & Yeager (2009) form: C¹⁰_N(U) = (a₀ + a₁ U + a₂ / U) × 10⁻³

Wind speed is clamped to U_min to avoid singularity in the a₂/U term.

surface_virtual_potential_temperature(
    T₀,
    p₀,
    constants,
    surface
) -> Any

Compute virtual potential temperature over a planar surface with surface temperature T₀ and surface pressure p₀,

\[θᵥ₀ = T₀ (1 + δᵛᵈ qᵛ⁺)\]

where $qᵛ⁺$ is the saturation specific humidity at the surface and $δᵛᵈ = Rᵛ/Rᵈ - 1$ (≈ 0.608 for water vapor in Earth's atmosphere; the actual value depends on the gas constants in constants).

dynamics_density_for_potential_temperature(
    dynamics::CompressibleDynamics
) -> Any

Return the density field for potential temperature diagnostics. For compressible dynamics, uses the actual density field.

dynamics_pressure_for_potential_temperature(
    dynamics::CompressibleDynamics
) -> Any

Return the pressure field for potential temperature diagnostics. For compressible dynamics, uses the actual pressure field.

dynamics_standard_pressure(
    dynamics::CompressibleDynamics
) -> Any

Return the standard pressure for potential temperature diagnostics.

buoyancy_forceᶜᶜᶜ(
    i,
    j,
    k,
    grid,
    dynamics::CompressibleDynamics,
    temperature,
    specific_moisture,
    microphysics,
    microphysical_fields,
    constants
) -> Any

Compute the buoyancy force density for compressible dynamics at cell center (i, j, k).

When a reference state is provided, the buoyancy force is computed as a perturbation:

\[ρ b = -g (ρ - ρ_r)\]

where $ρ_r$ is the reference density in discrete hydrostatic balance. This eliminates the $O(Δz^2)$ truncation error from the near-cancellation of $∂p/∂z$ and $gρ$, which is essential for stability with acoustic substepping at large time steps.

Without a reference state, the full gravitational force $-gρ$ is used.

compute_auxiliary_dynamics_variables!(
    model::AtmosphereModel{<:CompressibleDynamics}
)

Compute temperature and pressure jointly for CompressibleModel.

For compressible dynamics with potential temperature thermodynamics, temperature and pressure are coupled via the ideal gas law and the potential temperature definition:

\[θ = T (p₀/p)^κ \quad \text{and} \quad p = ρ R^m T\]

Eliminating the circular dependency gives the direct formula:

\[T = θ^γ \left(\frac{ρ R^m}{p₀}\right)^{γ-1}\]

where $γ = c_p / c_v$ is the heat capacity ratio. Once temperature is known, pressure is computed from the ideal gas law $p = ρ R^m T$.

This joint computation is necessary because, unlike anelastic dynamics where pressure comes from a reference state, compressible dynamics requires solving for both temperature and pressure simultaneously.

compute_dynamics_tendency!(
    model::AtmosphereModel{<:CompressibleDynamics}
)

Compute the density tendency for compressible dynamics using the continuity equation.

The density evolves according to:

\[\partial_t \rho = -\boldsymbol{\nabla \cdot} (\rho \boldsymbol{u})\]

Since momentum ρu is already available, this is simply the negative divergence of momentum.

compute_pressure_correction!(
    _::AtmosphereModel{<:CompressibleDynamics},
    Δt
)

No-op for CompressibleDynamics - pressure is computed diagnostically from the equation of state.

default_timestepper(
    dynamics::CompressibleDynamics
) -> Symbol

Return the default timestepper for CompressibleDynamics based on its time_discretization.

• SplitExplicitTimeDiscretization: Returns :AcousticSSPRungeKutta3 for acoustic substepping
• ExplicitTimeStepping: Returns :SSPRungeKutta3 for standard explicit time-stepping

dynamics_density(dynamics::CompressibleDynamics) -> Any

Return the prognostic density field for CompressibleDynamics.

dynamics_pressure(dynamics::CompressibleDynamics) -> Any

Return the pressure field for CompressibleDynamics. Pressure is computed diagnostically from the equation of state.

dynamics_pressure_solver(
    dynamics::CompressibleDynamics,
    grid
)

Return nothing for CompressibleDynamics - no pressure solver is needed. Pressure is computed directly from the equation of state.

dynamics_prognostic_fields(
    dynamics::CompressibleDynamics
) -> NamedTuple{(:ρ,), <:Tuple{Any}}

Return prognostic fields specific to compressible dynamics. Returns the density field as a prognostic variable.

make_pressure_correction!(
    _::AtmosphereModel{<:CompressibleDynamics},
    Δt
)

No-op for CompressibleDynamics - no pressure projection is needed.

materialize_dynamics(
    dynamics::CompressibleDynamics,
    grid,
    boundary_conditions,
    thermodynamic_constants
) -> Union{CompressibleDynamics{_A, D, P, _B, Nothing} where {_A, D<:(Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), P<:(Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), _B}, CompressibleDynamics{_A, D, P, _B, RS} where {_A, D<:(Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), P<:(Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), _B, RS<:(ExnerReferenceState{_A, FP, FD, FE} where {_A, FP<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), FD<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), FE<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B})})}}

Materialize a stub CompressibleDynamics into a full dynamics object with density and pressure fields.

mean_pressure(dynamics::CompressibleDynamics) -> Any

Return the mean (reference) pressure for CompressibleDynamics. For compressible dynamics, there is no separate mean pressure - returns the full pressure field.

pressure_anomaly(dynamics::CompressibleDynamics) -> Int64

Return the pressure anomaly for CompressibleDynamics. For compressible dynamics, there is no decomposition - returns zero.

standard_pressure(dynamics::CompressibleDynamics) -> Any

Return the standard pressure for potential temperature calculations.

surface_pressure(dynamics::CompressibleDynamics) -> Any

Return a standard surface pressure for boundary condition regularization. For compressible dynamics, uses the standard atmospheric pressure (101325 Pa).

total_pressure(dynamics::CompressibleDynamics) -> Any

Return the total pressure for CompressibleDynamics, in Pa.

compute_acoustic_substeps(
    grid,
    Δt,
    thermodynamic_constants
) -> Any

Compute the number of acoustic substeps from the horizontal acoustic CFL condition.

Uses a conservative sound speed estimate ℂᵃᶜ = √(γ Rᵈ Tᵣ) with Tᵣ = 300 K (giving ℂᵃᶜ ≈ 347 m/s) and the minimum horizontal grid spacing. The vertical CFL is not needed because the w-π' coupling is vertically implicit.

Following CM1, the substep count satisfies Δτ · ℂᵃᶜ / Δx_min ≤ 1 where Δτ = Δt / N is the acoustic substep size. A safety factor of 1.2 is applied to ensure stability with the forward-backward splitting.

convert_slow_tendencies!(substepper, model)

Convert slow momentum tendencies (Gˢρu, Gˢρv, Gˢρw) to slow velocity tendencies (uten, vten, wten) and slow pressure tendency (Gˢπ).

The velocity tendency is: uten ≈ Gˢρu / ρ The pressure tendency is: Gˢπ = -u · ∇π

These are frozen during the acoustic substep loop.

implicit_w_solve!(w, substepper, model, Δτ, π′_forcing)

Solve the vertically implicit w-π' system, then update w.

Instead of solving a tridiagonal system for w (which requires face-indexed arrays), we solve for π' at cell centers (matching the solver dimensions), then back-solve for w.

The approach:

1. Substitute $w⁺[k] = w[k] + Δτ Gˢw[k] - Δτ (cᵖ θᵥ / Δz) δz(π'⁺)$ into the pressure equation $π'⁺ = π' + π'_{forcing} - α Δτ S ∂w⁺/∂z$
2. This gives a tridiagonal system in π'⁺ at center locations
3. After solving for π'⁺, back-solve for w⁺ from the new pressure gradient

recover_full_fields!(model, substepper, U⁰, Δt_stage)

Recover full fields from Exner pressure acoustic variables.

After the acoustic loop, convert the updated velocity and Exner pressure back to Breeze's prognostic variables (ρ, ρu, ρv, ρw, ρθ).

For WS-RK3: The recovery uses U⁰ as the base state and adds the change computed by the acoustic loop. The velocity fields were modified in-place during the loop, so we need to compute the change and apply the WS-RK3 formula.

recover_full_fields_ssp!(model, substepper, α, U⁰, Δt)

SSP-RK3 recovery: $U_{new} = (1 - α) U⁰ + α U_{acoustic}$

Uses nonlinear recovery from π' to ρθ via the equation of state: $ρθ_{acoustic} = (p_{st}/R_d) (π_{ref} + π'_{final})^{c_v/R}$

This is exact for dry air and avoids linearization errors that could accumulate when the acoustic loop runs many substeps.

solve_for_pressure!(
    _::AtmosphereModel{<:CompressibleDynamics}
)

No-op for CompressibleDynamics - pressure is computed from the equation of state, not solved.

materialize_microphysical_fields(
    _::Breeze.Microphysics.DCMIP2016KesslerMicrophysics,
    grid,
    boundary_conditions
) -> NamedTuple{(:ρqᶜˡ, :ρqʳ, :qᵛ, :qᶜˡ, :qʳ, :precipitation_rate, :𝕎ʳ), <:Tuple{Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Nothing, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}}

Create and return the microphysical fields for the Kessler scheme.

Prognostic Fields (Density-Weighted)

• ρqᶜˡ: Density-weighted cloud liquid mass fraction.
• ρqʳ: Density-weighted rain mass fraction.

Diagnostic Fields (Mass Fractions)

• qᵛ: Water vapor mass fraction, diagnosed as $q^v = q^t - q^{cl} - q^r$.
• qᶜˡ: Cloud liquid mass fraction ($kg/kg$).
• qʳ: Rain mass fraction ($kg/kg$).
• precipitation_rate: Surface precipitation rate ($m/s$), defined as $q^r imes v^t_{rain}$ to match one-moment microphysics.
• 𝕎ʳ: Rain terminal velocity ($m/s$).

maybe_adjust_thermodynamic_state(
    𝒰,
    _::Breeze.Microphysics.DCMIP2016KesslerMicrophysics,
    qᵗ,
    constants
) -> Any

Return the thermodynamic state without adjustment.

The Kessler scheme performs its own saturation adjustment internally via the kernel.

microphysical_tendency(
    _::Breeze.Microphysics.DCMIP2016KesslerMicrophysics,
    name,
    ρ,
    ℳ,
    𝒰,
    constants
) -> Any

Return zero tendency.

All microphysical source/sink terms are applied directly to the prognostic fields via the microphysics_model_update! kernel, bypassing the standard tendency interface.

microphysical_velocities(
    _::Breeze.Microphysics.DCMIP2016KesslerMicrophysics,
    μ,
    name
)

Return nothing.

Rain sedimentation is handled internally by the kernel rather than through the advection interface.

microphysics_model_update!(
    microphysics::Breeze.Microphysics.DCMIP2016KesslerMicrophysics,
    model
)

Apply the Kessler microphysics to the model.

This function launches a kernel that processes each column independently, with rain sedimentation subcycling.

The kernel handles conversion between mass fractions and mixing ratios internally for efficiency. Water vapor is diagnosed from $q^v = q^t - q^{cl} - q^r$.

microphysics_model_update!(
    microphysics::Breeze.Microphysics.DCMIP2016KesslerMicrophysics,
    model::AtmosphereModel{<:ParcelDynamics}
)

Apply DCMIP2016 Kessler microphysics to a parcel model.

For a Lagrangian parcel, the microphysics processes are:

1. Autoconversion: Cloud water → rain when cloud exceeds threshold
2. Accretion: Rain + cloud → rain (collection)
3. Saturation adjustment: Vapor ↔ cloud to maintain equilibrium
4. Rain evaporation: Rain → vapor in subsaturated air

Note: Rain sedimentation is not applicable to a Lagrangian parcel since the parcel is a closed system (rain does not fall out of the parcel).

precipitation_rate(
    model,
    _::Breeze.Microphysics.DCMIP2016KesslerMicrophysics,
    _::Val{:liquid}
) -> Any

Return the liquid precipitation rate field for the DCMIP2016 Kessler microphysics scheme.

The precipitation rate is computed internally by the Kessler kernel and stored in μ.precipitation_rate. It is defined as $q^r imes v^t_{rain}$ (rain mass fraction times terminal velocity), matching the one-moment microphysics definition. Units are m/s.

This implements the Breeze precipitation_rate(model, phase) interface, allowing the DCMIP2016 Kessler scheme to integrate with Breeze's standard diagnostics.

prognostic_field_names(
    _::Breeze.Microphysics.DCMIP2016KesslerMicrophysics
) -> Tuple{Symbol, Symbol}

Return the names of prognostic microphysical fields for the Kessler scheme.

Fields

• :ρqᶜˡ: Density-weighted cloud liquid mass fraction ($kg/m^3$).
• :ρqʳ: Density-weighted rain mass fraction ($kg/m^3$).

surface_precipitation_flux(
    model,
    _::Breeze.Microphysics.DCMIP2016KesslerMicrophysics
) -> Field{LX, LY, Nothing, O, G, I, D, T, B, Oceananigans.Fields.FieldStatus{Float64}} where {LX, LY, O, G, I, D, T, B}

Return the surface precipitation flux field for the DCMIP2016 Kessler microphysics scheme.

The surface precipitation flux is $ ho q^r v^t_{rain}$ at the surface, matching the one-moment microphysics definition. Units are kg/m²/s.

This implements the Breeze surface_precipitation_flux(model) interface.

cloud_to_rain_production(rᶜˡ, rʳ, Δt, microphysics)

Compute cloud-to-rain production rate from autoconversion and accretion (Klemp & Wilhelmson 1978, eq. 2.13).

This implements the combined effect of:

• Autoconversion: Cloud water spontaneously converting to rain when rᶜˡ > rᶜˡ★
• Accretion: Rain collecting cloud water as it falls

The formula uses an implicit time integration for numerical stability.

condensation_rate(
    qᵛ,
    qᵛ⁺,
    qᶜˡ,
    T,
    ρ,
    q,
    τᶜˡ,
    constants
) -> Any

Compute the condensation/evaporation rate for cloud liquid water in a relaxation-to-saturation model.

This returns the rate of change of cloud liquid mass fraction (kg/kg/s). Positive values indicate condensation; negative values indicate evaporation. Evaporation is limited by the available cloud liquid.

deposition_rate(
    qᵛ,
    qᵛ⁺ⁱ,
    qᶜⁱ,
    T,
    ρ,
    q,
    τᶜⁱ,
    constants
) -> Any

Compute the deposition/sublimation rate for cloud ice in a relaxation-to-saturation model.

This returns the rate of change of cloud ice mass fraction (kg/kg/s). Positive values indicate deposition; negative values indicate sublimation. Sublimation is limited by the available cloud ice.

ice_thermodynamic_adjustment_factor(
    qᵛ⁺ⁱ,
    T,
    q,
    constants
) -> Any

Compute the thermodynamic adjustment factor Γ used in relaxation-to-saturation deposition/sublimation tendencies (ice analogue of thermodynamic_adjustment_factor).

mass_fractions_to_mixing_ratios(
    qᵗ,
    ρqᶜˡ,
    ρqʳ,
    ρ
) -> Tuple{Any, Any, Any}

Convert from mass fractions to mixing ratios.

Returns (rᵛ, rᶜˡ, rʳ) mixing ratios for use in Kessler physics.

mixing_ratios_to_mass_fractions(
    rᵛ,
    rᶜˡ,
    rʳ
) -> NTuple{4, Any}

Convert from mixing ratios back to mass fractions.

Returns (qᵛ, qᶜˡ, qʳ, qᵗ).

step_kessler_microphysics(
    rᵛ,
    rᶜˡ,
    rʳ,
    Δr𝕎,
    T,
    ρ,
    p,
    Δt,
    microphysics,
    constants,
    f₅,
    δT,
    FT
) -> NTuple{4, Any}

Apply one Kessler microphysics step: autoconversion, accretion, saturation adjustment, rain evaporation, and condensation.

Δr𝕎 is the sedimentation flux divergence (zero for parcel models).

Returns (rᵛ, rᶜˡ, rʳ, Δrˡ).

thermodynamic_adjustment_factor(qᵛ⁺, T, q, constants) -> Any

Compute the thermodynamic adjustment factor Γ used in relaxation-to-saturation condensation/evaporation tendencies.

compute_boussinesq_adjustment_temperature(
    𝒰₀::Breeze.Thermodynamics.LiquidIcePotentialTemperatureState{FT},
    constants
) -> Any

Return the temperature $T$ corresponding to thermodynamic equilibrium between the specific humidity and liquid mass fractions of the input thermodynamic state 𝒰₀, wherein the specific humidity is equal to or less than the saturation specific humidity at the given conditions and affiliated with theromdynamic constants constants.

The saturation equilibrium temperature satisfies the nonlinear relation

\[θ = [1 - ℒˡᵣ qˡ / (cᵖᵐ T)] T / Π ,\]

with $ℒˡᵣ$ the latent heat at the reference temperature $Tᵣ$, $cᵖᵐ$ the mixture specific heat, $Π$ the Exner function, $qˡ = \max(0, qᵗ - qᵛ⁺)$ the condensate specific humidity, $qᵗ$ is the total specific humidity, and $qᵛ⁺$ is the saturation specific humidity.

The saturation equilibrium temperature is thus obtained by solving $r(T) = 0$, where

\[r(T) ≡ T - θ Π - ℒˡᵣ qˡ / cᵖᵐ .\]

Solution of $r(T) = 0$ is found via the secant method.

initialize_parcel_state!(state, z₀, x₀, y₀, model)

Initialize the parcel state by interpolating environmental conditions at the given position.

Reconstruct a thermodynamic state with a new conserved variable value and updated z, p.

set_moisture_from_relative_humidity!(
    qᵗ_field,
    ℋ,
    T_field,
    ρ_field,
    constants
)

Set specific humidity field from relative humidity, computing qᵗ = ℋ * qᵛ⁺(T, ρ).

set_temperature_from_potential_temperature!(
    T_field,
    θ,
    p_field,
    pˢᵗ,
    constants
)

Set temperature field from potential temperature, using proper thermodynamic relations.

ssp_rk3_microphysics_substep(
    _::Nothing,
    _::Nothing,
    _::Nothing,
    Δt,
    α
)

Apply SSP RK3 substep formula to microphysics prognostic variables.

set!(
    model::AtmosphereModel{<:ParcelDynamics};
    T,
    θ,
    ρ,
    p,
    qᵗ,
    ℋ,
    u,
    v,
    w,
    x,
    y,
    z
)

Set the environmental profiles and initial parcel state for a ParcelModel.

Environmental profiles are set on the model's fields (temperature, density, pressure, velocities). The parcel is initialized at the specified position with environmental conditions interpolated at that height.

Keyword Arguments

Thermodynamic profiles (provide one of T or θ):

• T: Temperature profile T(z) [K] - function, array, Field, or constant
• θ: Potential temperature profile θ(z) [K] - function, array, or constant. If provided, T is computed from θ and p using thermodynamic relations.
• ρ: Density profile ρ(z) [kg/m³] - function, array, Field, or constant
• p: Pressure profile p(z) [Pa] - function, array, Field, or constant

Moisture (provide one of qᵗ or ℋ):

• qᵗ: Specific humidity profile qᵗ(z) [kg/kg] - function, array, or constant (default: 0)
• ℋ: Relative humidity profile ℋ(z) [0-1] - function, array, or constant. If provided, qᵗ is computed as qᵗ = ℋ * qᵛ⁺(T, ρ).

Velocities:

• u: Zonal velocity u(z) [m/s] - function, array, or constant (default: 0)
• v: Meridional velocity v(z) [m/s] - function, array, or constant (default: 0)
• w: Vertical velocity w(z) [m/s] - function, array, or constant (default: 0)

Parcel position:

• x: Initial parcel x-position [m] (default: 0)
• y: Initial parcel y-position [m] (default: 0)
• z: Initial parcel height [m] (required to initialize parcel state)

time_step!(
    model::AtmosphereModel{<:ParcelDynamics},
    Δt;
    callbacks
)

Advance the parcel model by one time step Δt using SSP RK3.

The SSP RK3 scheme Shu and Osher (1988) is:

\[u^{(1)} = u^{(0)} + Δt L(u^{(0)}) u^{(2)} = \frac{3}{4} u^{(0)} + \frac{1}{4} u^{(1)} + \frac{1}{4} Δt L(u^{(1)}) u^{(3)} = \frac{1}{3} u^{(0)} + \frac{2}{3} u^{(2)} + \frac{2}{3} Δt L(u^{(2)})\]

This scheme has CFL coefficient = 1 and is TVD (total variation diminishing).

update_state!(model::AtmosphereModel{<:ParcelDynamics}; ...)
update_state!(
    model::AtmosphereModel{<:ParcelDynamics},
    callbacks;
    compute_tendencies
)

Update the parcel model state, computing tendencies and auxiliary variables.

This function is called at the beginning of each time step and after each substep in multi-stage time steppers. It mirrors the role of update_state! for AtmosphereModel and consolidates all state-dependent computations:

1. Compute position tendencies (Gx, Gy, Gz) from environmental velocity profiles
2. Any other auxiliary state computations (currently none)

Keyword Arguments

• compute_tendencies: If true (default), compute tendencies for prognostic variables.

diagnose_thermodynamic_state(
    i,
    j,
    k,
    grid,
    formulation::LiquidIcePotentialTemperatureFormulation,
    dynamics,
    q
) -> Breeze.Thermodynamics.LiquidIcePotentialTemperatureState

Build a LiquidIcePotentialTemperatureState at grid point (i, j, k) from the given formulation, dynamics, and pre-computed moisture mass fractions q.

set_thermodynamic_variable!(
    model::AtmosphereModel{<:Any, <:LiquidIcePotentialTemperatureFormulation},
    _::Val{:T},
    value
)

Set the thermodynamic state from in-situ temperature $T$.

The temperature is converted to liquid-ice potential temperature θˡⁱ using the relation between $T$ and θˡⁱ` that accounts for the moisture distribution.

For unsaturated air (no condensate), this simplifies to $θ = T / Π$ where $Π$ is the Exner function.

static_energy_density(model::PotentialTemperatureModel)

Return the static energy density as a Field with boundary conditions that return energy fluxes when used with BoundaryConditionOperation.

For LiquidIcePotentialTemperatureFormulation, the prognostic variable is potential temperature density ρθ. This function converts the ρθ boundary conditions to energy flux boundary conditions by multiplying by the mixture heat capacity cᵖᵐ.

diagnose_thermodynamic_state(
    i,
    j,
    k,
    grid,
    formulation::StaticEnergyFormulation,
    dynamics,
    q
) -> Breeze.Thermodynamics.StaticEnergyState

Build a StaticEnergyState at grid point (i, j, k) from the given formulation, dynamics, and pre-computed moisture mass fractions q.

set_thermodynamic_variable!(
    model::AtmosphereModel{<:Any, <:StaticEnergyFormulation},
    _::Val{:T},
    value
)

Set the thermodynamic state from temperature $T$.

The temperature is converted to static energy $e$ using the relation:

\[e = cᵖᵐ T + g z - ℒˡ qˡ - ℒⁱ qⁱ .\]

enforce_discrete_hydrostatic_balance!(pᵣ, ρᵣ, grid, g)

Recompute the reference pressure pᵣ from the reference density ρᵣ by discrete upward integration, ensuring that the discrete hydrostatic balance

\[\frac{p_{ref}[k] - p_{ref}[k-1]}{Δz} + g \frac{ρ_{ref}[k] + ρ_{ref}[k-1]}{2} = 0\]

holds exactly at every interior z-face. This guarantees that reference-state subtraction in the pressure gradient and buoyancy cancels to machine precision, eliminating the $O(Δz^2)$ truncation error that would otherwise dominate the momentum tendency for nearly-hydrostatic flows.

numerically_integrated_hydrostatic_density(z, p₀, θ_func, pˢᵗ, constants)

Compute the dry hydrostatic density at height z from the numerically integrated pressure and the given potential temperature profile θ_func(z).

numerically_integrated_hydrostatic_pressure(z, p₀, θ_func, pˢᵗ, constants)

Compute the dry hydrostatic pressure at height z by numerically integrating dp/dz = -gρ from z=0, where ρ = p/(Rᵈ T) and T = θ(z) (p/pˢᵗ)^κ.

This function handles non-uniform potential temperature profiles θ(z) for which the closed-form adiabatic solution does not apply. Uses 1000 midpoint integration steps.

acoustic_rk3_substep!(model, Δt, β)

Apply a Wicker-Skamarock RK3 substep with acoustic substepping.

The acoustic substep loop handles momentum, density, and the thermodynamic variable (ρθ or ρe). The substep size is constant $Δτ = Δt/N$ across all stages, with the substep count varying as $Nτ = \mathrm{round}(β N)$. Remaining scalars (tracers) are updated with standard RK3.

acoustic_scalar_substep!(
    model,
    kernel!,
    Δt_implicit,
    kernel_args...
)

Update non-acoustic scalar fields (moisture, tracers) using the given kernel.

Iterates over prognostic fields, skipping the first 5 (ρ, ρu, ρv, ρw, ρθ) which are handled by the acoustic substep loop. For each remaining field, launches kernel! with the provided kernel_args and applies the implicit diffusion step.

acoustic_ssp_rk3_substep!(model, Δt, α, stage)

Apply an SSP RK3 substep with acoustic substepping.

The acoustic substep loop handles momentum, density, and the thermodynamic variable (ρθ or ρe). Remaining scalars (tracers) are updated using standard SSP RK3 with time-averaged velocities from the acoustic loop.

compute_slow_momentum_tendencies!(model)

Compute slow momentum tendencies (advection, Coriolis, turbulence, forcing).

The pressure gradient and buoyancy are excluded using SlowTendencyMode. These "fast" terms are handled by the acoustic substep loop, which resolves the acoustic CFL through substepping with constant $Δτ = Δt/N$.

compute_slow_scalar_tendencies!(model)

Compute slow tendencies for density and thermodynamic variable.

In the perturbation-variable approach, the slow tendencies are simply the full RHS $R^t$ evaluated at the stage-level state. No correction is needed because the acoustic loop advances perturbation variables, not full fields.

• $G^s_ρ = -\boldsymbol{∇·m}^t$: full density tendency (continuity equation)
• $G^s_χ$: full thermodynamic tendency (advection + physics)

time_step!(
    model::AtmosphereModel{<:CompressibleDynamics, <:Any, <:Any, <:AcousticRungeKutta3},
    Δt;
    callbacks
)

Step forward model one time step Δt with Wicker-Skamarock RK3 and acoustic substepping.

The algorithm follows Wicker and Skamarock (2002):

• Outer loop: 3-stage RK3 with stage fractions Δt/3, Δt/2, Δt
• Inner loop: Acoustic substeps for fast (pressure) tendencies

Each RK stage:

1. Compute slow tendencies (advection, Coriolis, diffusion only — PGF/buoyancy in acoustic loop)
2. Execute acoustic substep loop for momentum and density (full PGF + buoyancy)
3. Update scalars using standard RK update with time-averaged velocities

time_step!(
    model::AtmosphereModel{<:CompressibleDynamics, <:Any, <:Any, <:AcousticSSPRungeKutta3},
    Δt;
    callbacks
)

Step forward model one time step Δt with SSP RK3 and acoustic substepping.

The algorithm is the Wicker-Skamarock scheme:

• Outer loop: 3-stage SSP RK3 for slow tendencies
• Inner loop: Acoustic substeps for fast (pressure) tendencies

Each RK stage:

1. Compute slow tendencies (advection, Coriolis, diffusion)
2. Execute acoustic substep loop for momentum and density
3. Update scalars using standard RK update with time-averaged velocities

time_step!(
    model::Oceananigans.AbstractModel{<:SSPRungeKutta3},
    Δt;
    callbacks
) -> Vector

Step forward model one time step Δt with the SSP RK3 method.

The algorithm is:

u^(1) = u^(0) + Δt L(u^(0))
u^(2) = 3/4 u^(0) + 1/4 u^(1) + 1/4 Δt L(u^(1))
u^(3) = 1/3 u^(0) + 2/3 u^(2) + 2/3 Δt L(u^(2))

where L above is the right-hand-side, e.g., ∂u/∂t = L(u).

RadiativeTransferModel(
    grid::Oceananigans.Grids.AbstractGrid,
    ::AllSkyOptics,
    constants::ThermodynamicConstants;
    background_atmosphere,
    surface_temperature,
    coordinate,
    epoch,
    surface_emissivity,
    direct_surface_albedo,
    diffuse_surface_albedo,
    surface_albedo,
    solar_constant,
    liquid_effective_radius,
    ice_effective_radius,
    ice_roughness
)

Construct an all-sky (gas + cloud) full-spectrum RadiativeTransferModel for the given grid.

This constructor requires that NCDatasets is loadable in the user environment because RRTMGP loads lookup tables from netCDF via an extension.

Keyword Arguments

• background_atmosphere: Background atmospheric gas composition (default: BackgroundAtmosphere{FT}()).
• surface_temperature: Surface temperature in Kelvin (required).
• coordinate: Tuple of (longitude, latitude) in degrees. If nothing (default), extracted from grid coordinates.
• epoch: Optional epoch for computing time with floating-point clocks.
• surface_emissivity: Surface emissivity, 0-1 (default: 0.98). Scalar.
• surface_albedo: Surface albedo, 0-1. Can be scalar or 2D field. Alternatively, provide both direct_surface_albedo and diffuse_surface_albedo.
• direct_surface_albedo: Direct surface albedo, 0-1. Can be scalar or 2D field.
• diffuse_surface_albedo: Diffuse surface albedo, 0-1. Can be scalar or 2D field.
• solar_constant: Top-of-atmosphere solar flux in W/m² (default: 1361)
• liquid_effective_radius: Model for cloud liquid effective radius in μm (default: ConstantRadiusParticles(10.0))
• ice_effective_radius: Model for cloud ice effective radius in μm (default: ConstantRadiusParticles(30.0))
• ice_roughness: Ice crystal roughness for cloud optics (1=smooth, 2=medium, 3=rough; default: 2)

RadiativeTransferModel(
    grid::Oceananigans.Grids.AbstractGrid,
    ::ClearSkyOptics,
    constants::ThermodynamicConstants;
    background_atmosphere,
    surface_temperature,
    coordinate,
    epoch,
    surface_emissivity,
    direct_surface_albedo,
    diffuse_surface_albedo,
    surface_albedo,
    solar_constant
)

Construct a clear-sky (gas-only) full-spectrum RadiativeTransferModel for the given grid.

This constructor requires that NCDatasets is loadable in the user environment because RRTMGP loads lookup tables from netCDF via an extension.

Keyword Arguments

• background_atmosphere: Background atmospheric gas composition (default: BackgroundAtmosphere{FT}()).
• surface_temperature: Surface temperature in Kelvin (required).
• coordinate: Tuple of (longitude, latitude) in degrees. If nothing (default), extracted from grid coordinates.
• epoch: Optional epoch for computing time with floating-point clocks.
• surface_emissivity: Surface emissivity, 0-1 (default: 0.98). Scalar.
• surface_albedo: Surface albedo, 0-1. Can be scalar or 2D field. Alternatively, provide both direct_surface_albedo and diffuse_surface_albedo.
• direct_surface_albedo: Direct surface albedo, 0-1. Can be scalar or 2D field.
• diffuse_surface_albedo: Diffuse surface albedo, 0-1. Can be scalar or 2D field.
• solar_constant: Top-of-atmosphere solar flux in W/m² (default: 1361)

RadiativeTransferModel(
    grid::Oceananigans.Grids.AbstractGrid,
    ::GrayOptics,
    constants::ThermodynamicConstants;
    optical_thickness,
    surface_temperature,
    coordinate,
    epoch,
    surface_emissivity,
    direct_surface_albedo,
    diffuse_surface_albedo,
    surface_albedo,
    solar_constant
)

Construct a gray atmosphere radiative transfer model for the given grid.

Keyword Arguments

• optical_thickness: Optical thickness parameterization (default: GrayOpticalThicknessOGorman2008(FT)).
• surface_temperature: Surface temperature in Kelvin (required).
• coordinate: Tuple of (longitude, latitude) in degrees. If nothing (default), extracted from grid coordinates.
• epoch: Optional epoch for computing time with floating-point clocks.
• surface_emissivity: Surface emissivity, 0-1 (default: 0.98). Scalar.
• surface_albedo: Surface albedo, 0-1. Can be scalar or 2D field. Alternatively, provide both direct_surface_albedo and diffuse_surface_albedo.
• direct_surface_albedo: Direct surface albedo, 0-1. Can be scalar or 2D field.
• diffuse_surface_albedo: Diffuse surface albedo, 0-1. Can be scalar or 2D field.
• solar_constant: Top-of-atmosphere solar flux in W/m² (default: 1361)

RRTMGPParameters(constants::ThermodynamicConstants)

Construct RRTMGPParameters from Breeze's ThermodynamicConstants.

update_radiation!(
    rtm::RadiativeTransferModel{<:Any, <:Any, <:Any, <:Any, <:BackgroundAtmosphere, <:RRTMGP.AtmosphericStates.AtmosphericState{<:Any, <:Any, <:Any, <:Any, <:Any, <:RRTMGP.AtmosphericStates.CloudState}},
    model
)

Update the all-sky (gas + cloud) full-spectrum radiative fluxes from the current model state.

update_radiation!(
    rtm::RadiativeTransferModel{<:Any, <:Any, <:Any, <:Any, <:BackgroundAtmosphere},
    model
)

Update the clear-sky full-spectrum radiative fluxes from the current model state.

update_radiation!(
    rtm::RadiativeTransferModel{<:Any, <:Any, <:Any, <:Any, Nothing},
    model
)

Update the radiative fluxes from the current model state.

This function:

1. Updates the RRTMGP atmospheric state from model fields (T, p)
2. Computes the solar zenith angle from the model clock and grid location
3. Solves the longwave and shortwave RTE
4. Copies the fluxes to Oceananigans fields for output

Sign convention: positive flux = upward, negative flux = downward.

copy_fluxes_to_fields!(
    rtm::RadiativeTransferModel{<:Any, <:Any, <:Any, <:Any, Nothing},
    grid
)

Copy RRTMGP flux arrays to Oceananigans ZFaceFields.

Applies sign convention:

• positive = upward
• negative = downward.

For the non-scattering shortwave solver, only the direct beam flux is computed.

rrtmgp_context(arch::CPU) -> Any

Create an RRTMGP-compatible ClimaComms context from an Oceananigans architecture.

update_rrtmgp_state!(
    rrtmgp_state::RRTMGP.AtmosphericStates.GrayAtmosphericState,
    model,
    surface_temperature
)

Update the RRTMGP GrayAtmosphericState arrays from model fields.

Grid staggering: layers vs levels

RRTMGP requires atmospheric state at both "layers" (cell centers) and "levels" (cell faces). This matches the finite-volume staggering used in Oceananigans:

                        ┌─────────────────────────────────────────────────┐
    z_lev[Nz+1] ━━━━━━━ │  level Nz+1 (TOA):  p_lev, t_lev, z_lev         │ ← extrapolated
                        └─────────────────────────────────────────────────┘
                        ┌─────────────────────────────────────────────────┐
                        │  layer Nz:  T[Nz], p_lay[Nz] = pᵣ[Nz]           │ ← from model
                        └─────────────────────────────────────────────────┘
    z_lev[Nz]   ━━━━━━━   level Nz:   p_lev, t_lev, z_lev                   ← interpolated
                        ┌─────────────────────────────────────────────────┐
                        │  layer Nz-1                                     │
                        └─────────────────────────────────────────────────┘
                                            ⋮
                        ┌─────────────────────────────────────────────────┐
                        │  layer 2                                        │
                        └─────────────────────────────────────────────────┘
    z_lev[2]    ━━━━━━━   level 2:    p_lev, t_lev, z_lev                   ← interpolated
                        ┌─────────────────────────────────────────────────┐
                        │  layer 1:   T[1], p_lay[1] = pᵣ[1]              │ ← from model
                        └─────────────────────────────────────────────────┘
    z_lev[1]    ━━━━━━━   level 1 (surface, z=0):  p_lev = p₀, t_lev      │ ← from reference state
                        ══════════════════════════════════════════════════
                                        GROUND (t_sfc)

Why the model must provide level values

RRTMGP is a general-purpose radiative transfer solver that operates on columns of atmospheric data. It does not interpolate from layers to levels internally because:

1. Boundary conditions: The surface (level 1) and TOA (level Nz+1) require boundary values that only the atmospheric model knows. For pressure, we use the reference state's surface_pressure at z=0. For the top, we extrapolate using the adiabatic hydrostatic formula.

2. Physics-appropriate interpolation: Different quantities need different interpolation methods. Pressure uses geometric mean (log-linear interpolation) because it varies exponentially with height. Temperature uses arithmetic mean.

3. Model consistency: The pressure profile must be consistent with the atmospheric model's reference state. RRTMGP has no knowledge of the anelastic approximation or the reference potential temperature θ₀.

Physics notes

Temperature: We use the actual temperature field T from the model state. This is the temperature that matters for thermal emission and absorption.

Pressure: In the anelastic approximation, pressure perturbations are negligible compared to the hydrostatic reference pressure. We use reference_state.pressure at cell centers, computed via adiabatic_hydrostatic_pressure(z, p₀, θ₀).

RRTMGP array layout

• Layer arrays (nlay, ncol): values at cell centers, layer 1 at bottom
• Level arrays (nlev, ncol): values at cell faces, level 1 at surface (z=0)
• nlev = nlay + 1

update_solar_zenith_angle!(
    sw_solver,
    coordinate::Tuple,
    grid,
    datetime
)

Update the solar zenith angle in the shortwave solver from the model clock.

Uses the datetime from clock.time and the grid's location (latitude/longitude) to compute the cosine of the solar zenith angle via celestial mechanics.

Does not support anything but single-column grids for now.

AerosolActivation{AP, AD, FT}

Aerosol activation parameters for two-moment microphysics.

Aerosol activation is the physical process that creates cloud droplets from aerosol particles when air becomes supersaturated. This struct bundles the parameters needed to compute the activation source term for cloud droplet number concentration.

Fields

• activation_parameters: [AerosolActivationParameters] from CloudMicrophysics.jl
• aerosol_distribution: Aerosol size distribution (modes with number, size, hygroscopicity)
• nucleation_timescale: Nucleation timescale [s] for converting activation deficit to rate (default: 1s)

References

• Abdul-Razzak, H. and Ghan, S.J. (2000). A parameterization of aerosol activation:
  1. Multiple aerosol types. J. Geophys. Res., 105(D5), 6837-6844.

MixedPhaseOneMomentState{FT} <: AbstractMicrophysicalState{FT}

Microphysical state for mixed-phase one-moment bulk microphysics.

Contains the local mixing ratios for cloud liquid, cloud ice, rain, and snow. This state is used for both saturation adjustment and non-equilibrium cloud formation in mixed-phase simulations.

Fields

• qᶜˡ: Cloud liquid mixing ratio (kg/kg)
• qᶜⁱ: Cloud ice mixing ratio (kg/kg)
• qʳ: Rain mixing ratio (kg/kg)
• qˢ: Snow mixing ratio (kg/kg)

OneMomentCloudMicrophysics(FT = Oceananigans.defaults.FloatType;
                           cloud_formation = NonEquilibriumCloudFormation(nothing, nothing),
                           categories = one_moment_cloud_microphysics_categories(FT),
                           precipitation_boundary_condition = nothing)

Return a OneMomentCloudMicrophysics microphysics scheme for warm-rain and mixed-phase precipitation.

The one-moment scheme uses CloudMicrophysics.jl 1M processes:

• Condensation/evaporation of cloud liquid (relaxation toward saturation)
• Autoconversion of cloud liquid to rain
• Accretion of cloud liquid by rain
• Terminal velocity for rain sedimentation

By default, non-equilibrium cloud formation is used, where cloud liquid is a prognostic variable that evolves via condensation/evaporation tendencies following Morrison and Grabowski (2008) (see Appendix A). The prognostic variables are ρqᶜˡ (cloud liquid mass density) and ρqʳ (rain mass density).

For equilibrium (saturation adjustment) cloud formation, pass:

using Breeze.Microphysics
cloud_formation = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())

# output
SaturationAdjustment{WarmPhaseEquilibrium, Float64}(0.001, Inf, WarmPhaseEquilibrium())

Keyword arguments

• precipitation_boundary_condition: Controls whether precipitation passes through the bottom boundary.
  • nothing (default): Rain exits through the bottom (open boundary)
  • ImpenetrableBoundaryCondition(): Rain collects at the bottom (zero terminal velocity at surface)

See the CloudMicrophysics.jl documentation for details.

References

• Morrison, H. and Grabowski, W. W. (2008). A novel approach for representing ice microphysics in models: Description and tests using a kinematic framework. J. Atmos. Sci., 65, 1528–1548. https://doi.org/10.1175/2007JAS2491.1

TwoMomentCategories{W, AP, LV, RV, AA}

Parameters for two-moment (Seifert and Beheng, 2006) warm-rain microphysics.

Fields

• warm_processes: Seifert and Beheng (2006) parameters bundling autoconversion, accretion, self-collection, breakup, evaporation, number adjustment, and size distribution parameters
• air_properties: AirProperties for thermodynamic calculations
• cloud_liquid_fall_velocity: StokesRegimeVelType for cloud droplet terminal velocity
• rain_fall_velocity: SB2006VelType or Chen2022VelTypeRain for raindrop terminal velocity
• aerosol_activation: AerosolActivation parameters for cloud droplet nucleation (or nothing to disable)

References

• Seifert, A. and Beheng, K. D. (2006). A two-moment cloud microphysics parameterization for mixed-phase clouds. Part 1: Model description. Meteorol. Atmos. Phys., 92, 45-66. https://doi.org/10.1007/s00703-005-0112-4
• Abdul-Razzak, H. and Ghan, S.J. (2000). A parameterization of aerosol activation:
  1. Multiple aerosol types. J. Geophys. Res., 105(D5), 6837-6844.

TwoMomentCloudMicrophysics(FT = Oceananigans.defaults.FloatType;
                           cloud_formation = NonEquilibriumCloudFormation(nothing, nothing),
                           categories = two_moment_cloud_microphysics_categories(FT),
                           precipitation_boundary_condition = nothing)

Return a TwoMomentCloudMicrophysics microphysics scheme for warm-rain precipitation using the Seifert and Beheng (2006) two-moment parameterization.

The two-moment scheme tracks both mass and number concentration for cloud liquid and rain, using CloudMicrophysics.jl 2M processes:

• Aerosol activation: Creates cloud droplets when supersaturation develops (enabled by default)
• Condensation/evaporation of cloud liquid (relaxation toward saturation)
• Autoconversion of cloud liquid to rain (mass and number)
• Accretion of cloud liquid by rain (mass and number)
• Cloud liquid self-collection (number only)
• Rain self-collection and breakup (number only)
• Rain evaporation (mass and number)
• Number adjustment to maintain physical mean particle mass bounds
• Terminal velocities (number-weighted and mass-weighted)

Non-equilibrium cloud formation is used, where cloud liquid mass and number are prognostic variables that evolve via condensation/evaporation, aerosol activation, and microphysical tendencies.

The prognostic variables are:

• ρqᶜˡ: cloud liquid mass density [kg/m³]
• ρnᶜˡ: cloud liquid number density [1/m³]
• ρqʳ: rain mass density [kg/m³]
• ρnʳ: rain number density [1/m³]

Aerosol Activation

Aerosol activation is enabled by default and provides the physical source term for cloud droplet number concentration. Without activation, cloud droplets cannot form. The default aerosol population represents typical continental conditions (~100 cm⁻³).

To customize the aerosol population, pass a custom categories with different aerosol_activation:

# Marine aerosol (fewer, more hygroscopic particles)
marine_mode = CMAM.Mode_κ(0.08e-6, 1.8, 50e6, (1.0,), (1.0,), (0.058,), (1.0,))
marine_activation = AerosolActivation(
    AerosolActivationParameters(Float64),
    CMAM.AerosolDistribution((marine_mode,))
)
categories = two_moment_cloud_microphysics_categories(aerosol_activation = marine_activation)
microphysics = TwoMomentCloudMicrophysics(categories = categories)

Keyword arguments

• cloud_formation: Cloud formation scheme (default: NonEquilibriumCloudFormation)
• categories: TwoMomentCategories containing SB2006 and aerosol activation parameters
• precipitation_boundary_condition: Controls whether precipitation passes through the bottom boundary.
  • nothing (default): Rain exits through the bottom (open boundary)
  • ImpenetrableBoundaryCondition(): Rain collects at the bottom (zero terminal velocity at surface)

See the CloudMicrophysics.jl 2M documentation for details on the Seifert and Beheng (2006) scheme.

References

• Seifert, A. and Beheng, K. D. (2006). A two-moment cloud microphysics parameterization for mixed-phase clouds. Part 1: Model description. Meteorol. Atmos. Phys., 92, 45-66. https://doi.org/10.1007/s00703-005-0112-4
• Abdul-Razzak, H. and Ghan, S.J. (2000). A parameterization of aerosol activation:
  1. Multiple aerosol types. J. Geophys. Res., 105(D5), 6837-6844.

WarmPhaseOneMomentState{FT} <: AbstractMicrophysicalState{FT}

Microphysical state for warm-phase one-moment bulk microphysics.

Contains the local mixing ratios needed to compute tendencies for cloud liquid and rain. This state is used for both saturation adjustment and non-equilibrium cloud formation in warm-phase (liquid only) simulations.

Fields

• qᶜˡ: Cloud liquid mixing ratio (kg/kg)
• qʳ: Rain mixing ratio (kg/kg)

WarmPhaseTwoMomentState{FT, V} <: AbstractMicrophysicalState{FT}

Microphysical state for warm-phase two-moment bulk microphysics.

Contains the local mixing ratios and number concentrations needed to compute tendencies for cloud liquid and rain following the Seifert-Beheng 2006 scheme.

Fields

• qᶜˡ: Cloud liquid mixing ratio (kg/kg)
• nᶜˡ: Cloud liquid number per unit mass (1/kg)
• qʳ: Rain mixing ratio (kg/kg)
• nʳ: Rain number per unit mass (1/kg)
• nᵃ: Aerosol number per unit mass (1/kg)
• velocities: NamedTuple of velocity components (; u, v, w) [m/s]. The vertical velocity w is used for aerosol activation.

ZeroMomentCloudMicrophysics(FT = Oceananigans.defaults.FloatType;
                            cloud_formation = SaturationAdjustment(FT),
                            τ_precip = 1000,
                            qc_0 = 5e-4,
                            S_0 = 0)

Return a ZeroMomentCloudMicrophysics microphysics scheme for warm-rain precipitation.

The zero-moment scheme removes cloud liquid water above a threshold at a specified rate:

• τ_precip: precipitation timescale in seconds (default: 1000 s)

and either

• S_0: supersaturation threshold (default: 0)
• qc_0: cloud liquid water threshold for precipitation (default: 5×10⁻⁴ kg/kg)

For more information see the CloudMicrophysics.jl documentation.

ZeroMomentBulkMicrophysics

Type alias for BulkMicrophysics with CloudMicrophysics 0M precipitation scheme.

The 0M scheme instantly removes precipitable condensate above a threshold. Interface is identical to non-precipitating microphysics except that maybe_adjust_thermodynamic_state calls CloudMicrophysics remove_precipitation first.

precipitation_rate(
    model,
    microphysics::BulkMicrophysics{<:Any, <:CloudMicrophysics.Parameters.Parameters0M},
    _::Val{:liquid}
)

Return a Field representing the liquid precipitation rate (rain rate) in kg/kg/s.

For zero-moment microphysics, this is the rate at which cloud liquid water is removed by precipitation: -dqᵗ/dt from the remove_precipitation function.

surface_precipitation_flux(
    model,
    microphysics::BulkMicrophysics{<:Any, <:Breeze.Microphysics.FourCategories{<:CloudMicrophysics.Parameters.CloudLiquid, <:CloudMicrophysics.Parameters.CloudIce, <:CloudMicrophysics.Parameters.Rain, <:CloudMicrophysics.Parameters.Snow, <:CloudMicrophysics.Parameters.CollisionEff, <:CloudMicrophysics.Parameters.Blk1MVelType, <:CloudMicrophysics.Parameters.AirProperties}}
) -> Field{LX, LY, Nothing, O, G, I, D, T, B, Oceananigans.Fields.FieldStatus{Float64}} where {LX, LY, O, G, I, D, T, B}

Return a 2D Field representing the precipitation flux at the bottom boundary.

The surface precipitation flux is wʳ * ρqʳ at k=1 (bottom face), representing the rate at which rain mass leaves the domain through the bottom boundary.

Units: kg/m²/s (positive = downward, out of domain)

Note: The returned value is positive when rain is falling out of the domain (the terminal velocity wʳ is negative, and we flip the sign).

surface_precipitation_flux(
    model,
    microphysics::BulkMicrophysics{<:Any, <:BreezeCloudMicrophysicsExt.TwoMomentCategories{<:CloudMicrophysics.Parameters.SB2006, <:CloudMicrophysics.Parameters.AirProperties, <:CloudMicrophysics.Parameters.StokesRegimeVelType}}
) -> Field{LX, LY, Nothing, O, G, I, D, T, B, Oceananigans.Fields.FieldStatus{Float64}} where {LX, LY, O, G, I, D, T, B}

Return a 2D Field representing the precipitation flux at the bottom boundary.

The surface precipitation flux is wʳ * ρqʳ at k=1 (bottom face), representing the rate at which rain mass leaves the domain through the bottom boundary.

Units: kg/m²/s (positive = downward, out of domain)

Note: The returned value is positive when rain is falling out of the domain (the terminal velocity wʳ is negative, and we flip the sign).

aerosol_activated_fraction(aerosol_activation, aps, ρ, ℳ, 𝒰, constants)

Compute the fraction of aerosol that activates given current thermodynamic conditions. Uses the maximum supersaturation to determine which aerosol modes activate.

aerosol_activation_mass_tendency(aerosol_activation, aps, ρ, ℳ, 𝒰, constants)

Compute the cloud liquid mass tendency from aerosol activation.

When aerosol particles activate to form cloud droplets, the newly formed droplets have a finite initial size given by the activation radius. This function computes the corresponding mass source term for cloud liquid water.

The activation radius is derived from Köhler theory:

\[r_{act} = \frac{2A}{3 S}\]

where $A = 2σ/(ρ_w R_v T)$ is the curvature parameter and $S$ is the instantaneous supersaturation. See eq. 19 in Abdul-Razzak et al. (1998).

The mass tendency is then:

\[\frac{dq^{cl}}{dt}_{act} = \frac{dN^{cl}}{dt}_{act} \cdot \frac{4π}{3} r_{act}^3 \frac{ρ_w}{ρ}\]

The activation rate is controlled by the nucleation timescale τⁿᵘᶜ stored in the AerosolActivation parameters (default: 1s).

Returns

Mass tendency for cloud liquid [kg/kg/s]

default_aerosol_activation(FT = Float64; τⁿᵘᶜ = 1)

Create a default AerosolActivation representing a typical continental aerosol population.

The default distribution is a single mode with:

• Mean dry radius: 0.05 μm (50 nm)
• Geometric standard deviation: 2.0
• Number concentration: 100 cm⁻³ (100 × 10⁶ m⁻³)
• Hygroscopicity κ: 0.5 (typical for ammonium sulfate)

Keyword arguments

• τⁿᵘᶜ: Nucleation timescale [s] for converting activation deficit to rate (default: 1s). Controls how quickly the cloud droplet number relaxes toward the target activated number.

This provides sensible out-of-the-box behavior for two-moment microphysics. Users can customize the aerosol population by constructing their own AerosolActivation.

Example

# Use default aerosol
microphysics = TwoMomentCloudMicrophysics()

# Custom aerosol: marine (fewer, larger particles)
marine_mode = CMAM.Mode_κ(0.08e-6, 1.8, 50e6, (1.0,), (1.0,), (0.058,), (1.0,))
marine_aerosol = AerosolActivation(
    AerosolActivationParameters(Float64),
    CMAM.AerosolDistribution((marine_mode,)),
    1  # τⁿᵘᶜ = 1s
)
microphysics = TwoMomentCloudMicrophysics(aerosol_activation = marine_aerosol)

# Disable aerosol activation (not recommended)
microphysics = TwoMomentCloudMicrophysics(aerosol_activation = nothing)

diffusional_growth_factor(aps::AirProperties, T, constants)

Compute the thermodynamic factor $G$ that controls the rate of diffusional growth of cloud droplets and rain drops.

The $G$ factor combines the effects of thermal conductivity and vapor diffusivity on phase change. It appears in the Mason equation for droplet growth:

\[\frac{dm}{dt} = 4π r G 𝒮\]

where $𝒮$ is supersaturation and $r$ is droplet radius.

This is a translation of CloudMicrophysics.Common.G_func_liquid using Breeze's thermodynamics instead of Thermodynamics.jl.

See Eq. (13.28) by Pruppacher & Klett (2010).

References

• Pruppacher, H. R., Klett, J. D. (2010). Microphysics of clouds and precipitation. Springer Netherlands. 2nd Edition

max_supersaturation_breeze(aerosol_activation, aps, ρ, ℳ, 𝒰, constants)

Compute the maximum supersaturation using the Abdul-Razzak and Ghan (2000) parameterization.

This is a translation of CloudMicrophysics.AerosolActivation.max_supersaturation that uses Breeze's thermodynamics instead of Thermodynamics.jl.

Arguments

• aerosol_activation: AerosolActivation containing activation parameters and aerosol distribution
• aps: AirProperties (thermal conductivity, vapor diffusivity)
• ρ: Air density [kg/m³]
• ℳ: Microphysical state containing updraft velocity and number concentrations
• 𝒰: Thermodynamic state
• constants: Breeze ThermodynamicConstants

Returns

Maximum supersaturation (dimensionless, e.g., 0.01 = 1% supersaturation)

rain_evaporation(rain_params, vel, aps, q, qʳ, ρ, T, constants)

Compute the rain evaporation rate (dqʳ/dt, negative for evaporation).

This is a translation of CloudMicrophysics.Microphysics1M.evaporation_sublimation that uses Breeze's internal thermodynamics instead of Thermodynamics.jl.

Arguments

• rain_params: Rain microphysics parameters (pdf, mass, vent)
• vel: Terminal velocity parameters
• aps: Air properties (kinematic viscosity, vapor diffusivity, thermal conductivity)
• q: MoistureMassFractions containing vapor, liquid, and ice mass fractions
• qʳ: Rain specific humidity
• ρ: Air density
• T: Temperature
• constants: Breeze ThermodynamicConstants

Returns

Rate of change of rain specific humidity (negative = evaporation)

rain_evaporation_2m(sb, aps, q, qʳ, ρ, Nʳ, T, constants)

Compute the two-moment rain evaporation rate returning both number and mass tendencies.

This is a translation of CloudMicrophysics.Microphysics2M.rain_evaporation that uses Breeze's internal thermodynamics instead of Thermodynamics.jl.

Arguments

• sb: SB2006 parameters containing pdf_r and evap
• aps: Air properties (kinematic viscosity, vapor diffusivity, thermal conductivity)
• q: MoistureMassFractions containing vapor, liquid, and ice mass fractions
• qʳ: Rain specific humidity [kg/kg]
• ρ: Air density [kg/m³]
• Nʳ: Rain number concentration [1/m³]
• T: Temperature [K]
• constants: Breeze ThermodynamicConstants

Returns

Named tuple (; evap_rate_0, evap_rate_1) where:

• evap_rate_0: Rate of change of number concentration [1/(m³·s)], negative for evaporation
• evap_rate_1: Rate of change of mass mixing ratio [kg/kg/s], negative for evaporation

two_moment_cloud_microphysics_categories(FT = Oceananigans.defaults.FloatType;
                                         warm_processes = SB2006(FT),
                                         air_properties = AirProperties(FT),
                                         cloud_liquid_fall_velocity = StokesRegimeVelType(FT),
                                         rain_fall_velocity = SB2006VelType(FT),
                                         aerosol_activation = default_aerosol_activation(FT))

Construct TwoMomentCategories with default Seifert-Beheng 2006 parameters and aerosol activation.

Keyword arguments

• warm_processes: SB2006 parameters for warm-rain microphysics
• air_properties: Air properties for thermodynamic calculations
• cloud_liquid_fall_velocity: Terminal velocity parameters for cloud droplets (Stokes regime)
• rain_fall_velocity: Terminal velocity parameters for rain drops
• aerosol_activation: Aerosol activation parameters (default: continental aerosol). Set to nothing to disable activation (not recommended for physical simulations).