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:

\[\boldsymbol{\nabla} ⋅ (ρᵣ \boldsymbol{u}) = 0\]

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

Return AnelasticDynamics representing incompressible fluid dynamics expanded about reference_state.

total_energy_density_name

The key under which an energy input — a surface heat flux, a diabatic heating rate — is supplied to AtmosphereModel, :ρE, along with its specific alias :E for forcings.

$E$ denotes total energy (total_energy), so the key names the physical input without committing to the variable that carries it: a boundary_conditions or forcing entry keyed ρE is routed onto whichever thermodynamic density the formulation evolves (see thermodynamic_density_name), converted as that variable requires — divided by $cᵖᵐ Π$ for $ρθ$, whether it arrives as a flux or as a forcing, and passed through unconverted for $ρs$, which is itself an energy per unit mass.

$s$ names static energy specifically and $e$ is reserved for turbulent kinetic energy, so neither doubles as the energy key: ρs is a valid key only when static energy is the prognostic thermodynamic variable.

total_moisture_density_name

The key under which a water input — a surface evaporative flux, a moisture source — is supplied to AtmosphereModel, :ρqᵗ, along with its specific alias :qᵗ for forcings.

$qᵗ$ denotes total moisture, so the key names the physical input without committing to the variable that carries it: a boundary_conditions or forcing entry keyed ρqᵗ is routed onto whichever moisture density the microphysics evolves (see moisture_prognostic_name). Unlike an energy input, no conversion is involved: water added to the prognostic moisture is water added to $qᵗ$ under every scheme, so the routing is a pure re-key.

The prognostic moisture name is scheme-dependent — :ρqᵛ for non-equilibrium cloud formation, :ρqᵉ for saturation adjustment — and changes with the cloud_formation option of BulkMicrophysics, so keying a surface flux by that name ties a setup to one scheme. ρqᵗ is never itself prognostic and therefore works for any of them.

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.

abstract type AbstractNegativeMoistureCorrection

Abstract supertype for negative moisture correction schemes.

See fix_negative_moisture! for details.

abstract type AbstractNumberConcentrationCategories

Abstract supertype for microphysics categories that track number concentrations (e.g. two-moment schemes, aerosol-aware schemes).

Subtypes opt in to number concentration corrections (orphan zeroing and clamping) in the negative moisture correction. Schemes should extend correction_number_mass_pairs and correction_number_fields for their specific prognostic number fields.

abstract type AbstractSolarPosition

Abstract supertype for solar-position specifications passed to RadiativeTransferModel. Concrete subtypes determine how cos(θ_z) is computed on each radiation update:

• ApparentSolarPosition — real-Earth time-varying, computed from the model clock and grid (or explicit) longitude/latitude.
• DiurnalSolarPosition — idealized diurnal cycle at a fixed latitude and declination, no calendar dependence.
• FixedCosineZenith — constant cos(θ_z), clock-independent.

struct AdiabaticBalancer{T, S}

Configuration for adiabatic (FV3 na_init) initialization, applied with balance_adiabatically!(model, balancer) or set!(model; balancer = AdiabaticBalancer(...)). Works for both CompressibleDynamics and AnelasticDynamics.

Keyword arguments

• time_stepping: the time discretization used for the balance excursion (the sponge is always stripped — it is irreversible). CompressibleDynamics only; ignored for AnelasticDynamics (which has a single projection-based scheme). Options:
  • default (DefaultTimeStepping()) — fully-explicit stepping. Memory-minimal (no acoustic substepper; only the aliased Gⁿ/U⁰ tendency storage) and cleanly reversible, but Δt is bounded by the vertical acoustic CFL.
  • nothing — reuse the model's native scheme (e.g. split-explicit), at the cost of rebuilding the acoustic substepper's scratch fields.
  • any time-discretization object — swapped in as-is.
• Δt: forward/backward step size. nothing (default) auto-derives the vertical-acoustic-CFL step acoustic_cfl_safety · Δz_min / c from the grid and analysis sound speed; pass a number to override.
• cycles: number of balance cycles (default 1).
• weight: nudging weight toward the analysis snapshot (default 2 → ⅓ dynamics + ⅔ analysis).
• with_moisture: if true (default) the moisture density ρqᵉ relaxes with the other prognostics. If false, ρqᵉ is snapshotted before the balance and restored after, so it is preserved exactly — reproducing a graft that returns only (ρ, ρu, ρv, ρw, ρθ).

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.

struct ApparentSolarPosition{C, E} <: AbstractSolarPosition

Time-varying apparent solar position. The cosine of the solar zenith angle is recomputed on each radiation update from the model clock and either the grid's $(λ, φ)$ coordinates (when coordinate === nothing, the default) or an explicit (longitude, latitude) tuple stored in coordinate.

When the model clock holds a floating-point time (in seconds), epoch::DateTime provides the absolute reference against which clock.time is resolved. With a DateTime clock, epoch is ignored.

Fields

• coordinate::Any: Observer longitude/latitude. Either nothing (use grid coordinates) or a (longitude, latitude) tuple in degrees.

• epoch::Any: DateTime anchor for floating-point clocks. Either nothing (requires a DateTime clock) or a DateTime.

ApparentSolarPosition(
;
    coordinate,
    epoch
) -> ApparentSolarPosition{Nothing, Nothing}

Construct an ApparentSolarPosition with optional coordinate and epoch.

julia> using Breeze, Datesjulia> ApparentSolarPosition()ApparentSolarPosition(coordinate=<from grid>, epoch=<from clock>)julia> ApparentSolarPosition(coordinate = (-70.9, 42.5))ApparentSolarPosition(coordinate=(-70.9, 42.5), epoch=<from clock>)julia> ApparentSolarPosition(epoch = DateTime(2024, 1, 1))ApparentSolarPosition(coordinate=<from grid>, epoch=2024-01-01T00:00:00)

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,
    particles
) -> AtmosphereModel{Dyn, Frm, Arc, Tst, Grd, Clk, Thm, Mom, Moi, Nothing, Tmp, Sol, Vel, Trc, Adv, Nothing, Frc, Nothing, Cnd, Nothing, Nothing, Nothing, Nothing} where {Dyn, Frm, Arc, Tst, Grd, Clk, Thm, Mom, Moi, Tmp, Sol, Vel, Trc, Adv, Frc, Cnd}

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.

• particles are Lagrangian particles to be advected with the flow, constructed with Oceananigans.LagrangianParticles. Particles are advected with the Cartesian velocities model.velocities once per time step, over the full Δt. Default: nothing. See the "Lagrangian particles" section of the documentation for details, including the treatment on terrain-following grids.

Example

julia> using Breezejulia> 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{N2, O2, CO2, CH4, N2O, CO, NO2, O3, CFC11, CFC12, CFC22, CCL4, CF4, HFC125, HFC134A, HFC143A, HFC23, HFC32}

Volume mixing ratios (VMR) for radiatively active gases. All values are dimensionless molar fractions.

RRTMGP supports spatially-varying VMR only for H₂O (computed from model moisture) and O₃. All other gases use global mean values.

Fields

• Constant gases (global mean only): N₂, O₂, CO₂, CH₄, N₂O, CO, NO₂
• Halocarbons: CFC₁₁, CFC₁₂, CFC₂₂, CCl₄, CF₄
• Hydrofluorocarbons: HFC₁₂₅, HFC₁₃₄ₐ, HFC₁₄₃ₐ, HFC₂₃, HFC₃₂
• Spatially-varying: O₃ - can be a constant or a function for height-dependent profiles

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.

The BackgroundAtmosphere constructor does not require a grid. When passed to RadiativeTransferModel, the O₃ field is materialized using the grid. This allows users to seamlessly switch between constant and function-based concentrations.

BackgroundAtmosphere(
;
    N₂,
    O₂,
    CO₂,
    CH₄,
    N₂O,
    CO,
    NO₂,
    O₃,
    CFC₁₁,
    CFC₁₂,
    CFC₂₂,
    CCl₄,
    CF₄,
    HFC₁₂₅,
    HFC₁₃₄ₐ,
    HFC₁₄₃ₐ,
    HFC₂₃,
    HFC₃₂
) -> BackgroundAtmosphere{Float64, Float64, Float64, Float64, Float64, Float64, Float64, typeof(standard_ozone_profile), Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64}

Construct a BackgroundAtmosphere with volume mixing ratios for radiatively active gases. All values are dimensionless molar fractions.

RRTMGP supports spatially-varying VMR only for H₂O and O₃. Other gases use global means.

• Constant gases: Specify as numbers
• O₃: Can be a Number or Function for height-dependent profiles

Keyword Arguments

• Constant gases: N₂, O₂, CO₂, CH₄, N₂O, CO, NO₂
• Halocarbons: CFC₁₁, CFC₁₂, CFC₂₂, CCl₄, CF₄
• Hydrofluorocarbons: HFC₁₂₅, HFC₁₃₄ₐ, HFC₁₄₃ₐ, HFC₂₃, HFC₃₂
• Spatially-varying: O₃ (can be Number or Function)

Defaults are approximate modern atmospheric values; halocarbons default to zero, and ozone defaults to standard_ozone_profile (pass O₃ = 0 for an ozone-free atmosphere). Note: H₂O is computed from the model's prognostic moisture field.

Example

julia> using Breezejulia> background = BackgroundAtmosphere(CO₂ = 400e-6)BackgroundAtmosphere with 6 active gases:  N₂ = 0.78084  O₂ = 0.20946  CO₂ = 400.0 ppm  CH₄ = 1.8 ppm  N₂O = 330.0 ppb  O₃ = standard_ozone_profile (generic function with 1 method)julia> tropical_ozone(z) = 30e-9 * (1 + z / 10000);julia> background = BackgroundAtmosphere(CO₂ = 400e-6, O₃ = tropical_ozone)BackgroundAtmosphere with 6 active gases:  N₂ = 0.78084  O₂ = 0.20946  CO₂ = 400.0 ppm  CH₄ = 1.8 ppm  N₂O = 330.0 ppb  O₃ = tropical_ozone (generic function with 1 method)

A callable that returns the advective timescale of a model restricted to the directions of formulation: HorizontalFormulation() counts only the horizontal advective CFL (dropping the vertical term), ThreeDimensionalFormulation() counts all three directions. Pass it to the cell_advection_timescale keyword of TimeStepWizard / conjure_time_step_wizard!, or as the timescale argument of CFL (CFL(Δt, CellAdvectionTimescale(...))), to control or monitor which directions bind the time step.

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 meters.

struct DefaultTemperatureSolver

Sentinel indicating that a formulation's temperature solver should be chosen by the dynamics: materialize_formulation replaces it with default_temperature_solver(dynamics).

DiurnalSolarPosition(; ...)
DiurnalSolarPosition(
    FT::DataType;
    latitude,
    declination,
    day_length,
    noon_offset
)

Construct a DiurnalSolarPosition with sensible defaults: perpetual equinox (declination = 0), 24-hour day (day_length = 86400 s), and noon at the start of the simulation (noon_offset = 0).

The positional argument FT controls the precision of the stored fields and defaults to Oceananigans.defaults.FloatType. Pass FT = Float32 (or set Oceananigans.defaults.FloatType = Float32) to run in Float32:

DiurnalSolarPosition(Float32, latitude = 30)

struct DiurnalSolarPosition{FT} <: AbstractSolarPosition

Idealized diurnal cycle with no annual variation and no calendar dependence. cos(θ_z) is computed analytically on each radiation update from the model clock (which must be numeric — seconds since the start of the run) as

\[\cos(θ_z) = \sin(φ) \sin(δ) + \cos(φ) \cos(δ) \cos(ω), \qquad ω = \frac{2π}{T_d} (t - t_{\text{noon}})\]

where $φ$ is the (fixed) observer latitude, $δ$ is the (fixed) solar declination, $T_d$ is the day length, and $t_{\text{noon}}$ is the simulation time at which local noon occurs. $ω = 0$ at noon and $ω = ±π$ at local midnight. The result is clamped to be non-negative.

Fields

• latitude::Any: Observer latitude (degrees).

• declination::Any: Solar declination (degrees). Zero is perpetual equinox; ±23.5 is perpetual solstice.

• day_length::Any: Rotation period (seconds). Default 86400 is the Earth day.

• noon_offset::Any: Simulation time (seconds) at which local noon occurs. Default 0.

Examples

Perpetual equinox at 30°N (default: 24-hour day, noon at $t = 0$):

julia> using Breezejulia> DiurnalSolarPosition(latitude = 30)DiurnalSolarPosition(latitude = 30.0°, declination = 0.0°, day_length = 86400.0 s, noon_offset = 0.0 s)

Perpetual June solstice at 45°N:

julia> using Breezejulia> DiurnalSolarPosition(latitude = 45, declination = 23.5)DiurnalSolarPosition(latitude = 45.0°, declination = 23.5°, day_length = 86400.0 s, noon_offset = 0.0 s)

Fast rotator with a 10-hour day, sun overhead, starting at sunrise:

julia> using Breezejulia> DiurnalSolarPosition(latitude = 0, day_length = 10 * 3600, noon_offset = 5 * 3600)DiurnalSolarPosition(latitude = 0.0°, declination = 0.0°, day_length = 36000.0 s, noon_offset = 18000.0 s)

struct FixedCosineZenith{FT} <: AbstractSolarPosition

Constant cosine of the solar zenith angle. The model clock has no effect on the sun position; the shortwave path length is fixed at $1 / \cos(θ_z)$ and the top-of-atmosphere downward shortwave flux is solar_constant * cos_zenith.

This is the appropriate choice for idealized studies (radiative-convective equilibrium, RCE intercomparisons) where a diurnal or annual mean is desired. Common values: $\cos(θ_z) = 0.5$ for diurnal mean at mid-latitudes, $\cos(θ_z) ≈ 0.41$ for the global annual mean.

Fields

• cos_zenith::Any: Cosine of the solar zenith angle. Should satisfy $0 ≤ \cos(θ_z) ≤ 1$ for the sun above the horizon.

Example

julia> using Breezejulia> FixedCosineZenith(0.5)FixedCosineZenith(cos_zenith = 0.5)

struct GrayOptics <: Breeze.AtmosphereModels.AbstractOptics

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

References

• O'Gorman, P. A. and Schneider, T. (2008). The hydrological cycle over a wide range of climates simulated with an idealized GCM. Journal of Climate, 21, 3815–3832.

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.

HydrostaticallyBalancedDensity(; surface_pressure = nothing)

Marker passed as the ρ value to set! to set the density in discrete moist hydrostatic balance with the just-set θˡⁱ/qᵛ, by per-column integration of the hydrostatic equation upward from the pressure at the bottom face of each column. For CompressibleDynamics.

When the dynamics carries an ExnerReferenceState, the default anchor is taken from it, so the balanced state and the reference it will be differenced against use the same per-column pressure. Without a reference, the anchor is obtained by reducing the dynamics' $z = 0$ datum to each column's bottom face along its current near-surface thermodynamic state. This matters on a terrain-following grid, where anchoring every column at one scalar instead would leave the cold start with no surface pressure gradient across the terrain.

surface_pressure overrides that anchor with a scalar applied to every column. On a terrain-following grid, prefer the default: a scalar cannot represent the terrain-following surface pressure, and one that disagrees with the reference reintroduces the inconsistency.

Unlike supplying a density field, this guarantees the initial column satisfies the discrete hydrostatic balance (pᵏ − pᵏ⁻¹)/Δz + g(ρᵏ + ρᵏ⁻¹)/2 = 0, so the cold start carries no spurious vertical pressure-gradient force. Combine with compute_reference_state = true (perturbation-form base state) and balancer (nonhydrostatic ρw spin-up) for a full one-call initialization.

The current column solve supports liquid-ice potential-temperature thermodynamics and vapor-only moisture. It rejects nonzero liquid/ice condensate because condensate heat capacity and latent corrections are not yet included in the column integration.

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.

Solar position

The solar_position keyword controls how the cosine of the solar zenith angle is obtained on each radiation update. See AbstractSolarPosition and its subtypes:

• ApparentSolarPosition (default) — time-varying, computed from the model clock and grid (or explicit) longitude/latitude. Supports DateTime clocks and floating-point clocks resolved against an epoch.
• FixedCosineZenith — constant cos(θ_z), clock-independent. Appropriate for idealized radiative-convective equilibrium studies.

Example

julia> using Breeze, Oceananigans.Units, RRTMGP, NCDatasetsjulia> 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⁻²├── solar_position: ApparentSolarPosition(coordinate=(0.0, 45.0), epoch=<from clock>)├── 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, GrayOptics(), ThermodynamicConstants();                              surface_temperature = 300,                              surface_albedo = 0.1,                              solar_position = FixedCosineZenith(0.5))RadiativeTransferModel├── solar_constant: 1361.0 W m⁻²├── solar_position: FixedCosineZenith(cos_zenith = 0.5)├── 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⁻²├── solar_position: ApparentSolarPosition(coordinate=(0.0, 45.0), epoch=<from clock>)├── surface_temperature: ConstantField(300.0) K├── surface_emissivity: ConstantField(0.98)├── direct_surface_albedo: ConstantField(0.1)└── diffuse_surface_albedo: ConstantField(0.1)

References

• O'Gorman, P. A. and Schneider, T. (2008). The hydrological cycle over a wide range of climates simulated with an idealized GCM. Journal of Climate, 21, 3815–3832.

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.

struct SpeciesBorrowing{VB} <: Breeze.AtmosphereModels.AbstractNegativeMoistureCorrection

Correct negative moisture produced by advection via same-level species borrowing.

At each grid cell, negative hydrometeors borrow from lighter species in the chain (e.g. rain <- cloud liquid <- vapor). Vertical redistribution of any remaining negative vapor is performed when vertical_borrowing is set to VerticalBorrowing.

For microphysics with number concentrations (categories that subtype AbstractNumberConcentrationCategories), orphaned number concentrations are zeroed and negative number concentrations are clamped after mass borrowing.

See fix_negative_moisture! for details.

Fields

• vertical_borrowing: nothing (default) or VerticalBorrowing() to enable vertical redistribution

struct VerticalBorrowing <: Breeze.AtmosphereModels.AbstractNegativeMoistureCorrection

Redistribute negative vapor vertically within each column via a top-to-bottom sweep that pushes deficits downward, followed by one bottom-to-top borrowing step if the bottom level is still negative.

This scheme can be used on its own to correct the moisture prognostic, or as the second phase of SpeciesBorrowing to clean up any vapor deficit that remains after same-level species borrowing.

Column-integrated moisture is conserved ($Δz$-weighted).

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]

advecting_momentum(
    model
) -> NamedTuple{(:ρu, :ρv, :ρw), <:Tuple{Any, Any, Any}}

Return the momentum tuple used for momentum advection transport and the continuity equation divergence.

For standard (non-terrain) models, this is model.momentum. For terrain-following coordinates, the vertical component ρw is replaced by the contravariant vertical momentum $\rho \tilde{w}$.

aerosol_field_names(
    microphysics
) -> Union{Tuple{}, Tuple{Symbol}}

Return the names of prognostic fields that carry aerosol populations for microphysics.

Schemes without prognostic aerosol return an empty tuple by default. Microphysics schemes with prognostic aerosol extend this interface so model components can retain or process those fields without depending on scheme-specific names.

balance_adiabatically!(
    model::AtmosphereModel,
    balancer::AdiabaticBalancer
) -> AtmosphereModel

Run adiabatic (FV3 na_init) initialization on model in place: spin the nonhydrostatic state (ρw and the pressure balance) into balance with the analysis fields. balancer is an AdiabaticBalancer (or true for the defaults / false for a no-op). Builds a stripped, memory-sharing twin via adiabatic_balance_twin and runs the low-level balance_adiabatically!(model; Δt, cycles, weight) on it, so the balanced state lands directly in model — no graft, no second field set.

balance_adiabatically!(
    model::AtmosphereModel;
    Δt,
    cycles,
    weight
)

Spin up a balanced vertical momentum ρw (and the nonhydrostatic pressure balance) consistent with model's initial (analysis) state, via FV3 adiabatic initialization (na_init).

Analyses (ERA5, GFS, …) supply the density, momentum, and thermodynamic state but cold-start the vertical velocity w at zero (hydrostatic), so the nonhydrostatic state is out of balance with the rest. Each of cycles cycles entails two symmetric forward/backward dynamics excursions at the same Δt. After each excursion — which lets ρw develop — the initial fields (every prognostic except ρw) are nudged back toward their t = 0 snapshot by the weighted mean

x ← (x + weight·x₀) / (1 + weight)

(default weight = 2 → ⅓ dynamics + ⅔ snapshot). ρw is never snapshotted or nudged, so the balance the excursion imprints on it is exactly what is kept. update_state! after each nudge rebuilds the diagnostics; the clock is reset to t = 0 on exit.

balance_adiabatically! performs adiabatic dynamics only. The caller must pass a model built without physics (microphysics = nothing), without an upper sponge, and without forcing — these run inside update_state!/time_step! and would corrupt the spin-up. Boundary conditions are not modified; pass a model whose boundaries are time-invariant so the symmetric excursion stays nearly reversible. The two-argument balance_adiabatically!(model, balancer) constructs such a model automatically.

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

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

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_microphysical_tendencies!(model) -> Any

Add microphysics tendency contributions to the model's Gⁿ fields.

This is the only entry point through which compute_tendencies! adds microphysical sources to the model's Gⁿ fields. Concrete implementations add methods on the two-argument helper compute_microphysical_tendencies!(microphysics, model).

The default implementation launches a single fused kernel that builds the microphysical state ℳ and thermodynamic state 𝒰 once per cell, then +=s the result of microphysical_tendency for each prognostic name into the corresponding G field. Schemes whose tendencies factor naturally per-name only need to extend microphysical_tendency.

Schemes whose tendencies bundle many process rates feeding multiple prognostics (e.g. mixed-phase non-equilibrium 1M, where ~14 process rates feed 5 prognostic tendencies) override this method directly to compute the bundle once per cell.

compute_pressure_correction!(model, Δt)

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

default_temperature_solver(dynamics)

Return the default solver for a formulation's temperature inversion given dynamics.

The need for an iterative inversion is dictated by the intersection of the dynamics and the thermodynamic formulation: the fallback returns nothing (closed-form, no iteration), and dynamics whose prognostic closure makes the inversion implicit (e.g. CompressibleDynamics with LiquidIcePotentialTemperatureFormulation, where temperature solves T = (ρRᵐT/pˢᵗ)^κ θ + ΔL/cᵖᵐ) extend this function to return an iterative solver.

dynamics_density(dynamics)

Return the coupling density — the density weighting the momentum (ρu = ρᵈ u) and the thermodynamic flux variable (ρθ = ρᵈ θ), and the divisor for diagnosing velocity (u = ρu/ρᵈ) and potential temperature (θ = ρθ/ρᵈ). It is the prognostic mass variable advanced by continuity.

• AnelasticDynamics: the time-independent reference density ρᵣ.
• CompressibleDynamics: the prognostic dry-air density ρᵈ.

The total air density ρ = ρᵈ + Σ ρˣ (dry air plus every water species) is a separate, diagnosed quantity — see total_density — used wherever total mass enters the physics: the moisture mass-fraction recovery (qˣ = ρˣ/ρ, so the thermodynamics stays in mass fractions), scalar and water advection, the equation of state, and buoyancy. The water densities (ρqᵛ, ρqˡ, …) are stored as partial densities (mass per volume), not coupling-weighted. On the anelastic core the two densities coincide (total_density === dynamics_density).

dynamics_pressure(dynamics)

Return the pressure field appropriate to the dynamical formulation, in Pa — the pressure entering the equation of state, buoyancy, and the thermodynamic tendencies.

For anelastic dynamics, this is the time-independent hydrostatic reference pressure $pᵣ(z)$, excluding the non-hydrostatic pressure anomaly that enforces the divergence constraint but does not perturb the thermodynamic state. For compressible dynamics, this is the diagnosed equation-of-state pressure. The anomaly and total-pressure counterparts are pressure_anomaly and total_pressure.

This is the pressure every physics parameterization should read, including radiation. Call the accessor rather than reaching into dynamics.reference_state directly: which of the two the reference state is varies by formulation. On the anelastic core it is the thermodynamic pressure, and dynamics_pressure returns it. On the compressible core (flat or terrain) it is a pressure-gradient device built once from a fixed profile, it does not track the thermodynamic state, and several dynamics types carry none at all. total_density is the density counterpart.

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.

initial_aerosol_number(microphysics) -> Any

Return the aerosol population stored in a microphysics scheme's native units.

The units are the scheme's own: a volumetric distribution returns [m⁻³], while a distribution specified per unit mass of air returns [kg⁻¹]. Use initial_aerosol_number_density to obtain the value a prognostic ρnᵃ holds, whichever basis a scheme uses. A scheme may report a population here and carry no ρnᵃ.

Returns 0 by default.

initial_aerosol_number_density(microphysics, ρ) -> Any

Return the default aerosol number density $ρ nᵃ$ [m⁻³] for a microphysics scheme, given the air density ρ (a field for grid models, a number for parcels).

This is the value set! writes into the prognostic field ρnᵃ when the scheme has one and the user supplies neither nᵃ nor ρnᵃ. It is derived from the aerosol size distribution stored in the microphysics scheme, so it stays consistent with the activation parameters. A scheme may report a population here and still have no field to write it to; aerosol_field_names(microphysics) == () is what says so.

Each scheme is responsible for the units of its own aerosol distribution: the density argument is here so that a scheme whose distribution is specified per unit mass [kg⁻¹], as PredictedParticlePropertiesMicrophysics is through AerosolMode.number_mixing_ratio, can return the $ρ$-weighted value that ρnᵃ holds, while a scheme whose distribution is already a volumetric concentration [m⁻³] ignores ρ. Breeze's convention throughout is that nᵃ = ρnᵃ / ρ is per unit mass [kg⁻¹]; a scheme that omits this scaling diagnoses nᵃ with a spurious inverse-density dependence.

By default, forwards to initial_aerosol_number, which returns 0 unless a scheme overrides it.

is_density_tendency_forcing(_) -> Bool

Return true if forcing produces a density-weighted tendency F_{ρϕ} directly (i.e., already includes the multiplication by ρ).

Forcings that return density tendencies must be supplied under their density-weighted key (e.g., ρθ, ρu) rather than the corresponding specific key (θ, u), because the specific-key dispatch wraps user values in SpecificForcing, which would multiply by ρ a second time. This trait is used by atmosphere_model_forcing to reject such misuses with a clear error.

Defaults to false. Extended for SubsidenceForcing and GeostrophicForcing in the Forcings module.

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, thermodynamic_constants)

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

Boundary conditions supplied under the energy key ρE (see total_energy_density_name) are routed onto the prognostic thermodynamic variable of formulation: for :LiquidIcePotentialTemperature they become ρθ boundary conditions, with flux BCs wrapped in EnergyFluxBoundaryCondition to divide by the local mixture heat capacity; for :StaticEnergy they pass through onto ρs unconverted.

Conditions supplied under the moisture key ρqᵗ (see total_moisture_density_name) are likewise routed onto the moisture density that microphysics evolves — ρqᵛ or ρqᵉ, depending on the scheme — without conversion, since water entering the prognostic moisture is water entering $qᵗ$ under any of them.

The dynamics argument provides the standard pressure that boundary conditions need at materialization time.

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

Nothing about the model state is captured here: surface fluxes read the pressure, density, temperature and moisture they need from the model field tuple at evaluation time.

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.

materialize_background_atmosphere(
    atm::BackgroundAtmosphere,
    grid
) -> BackgroundAtmosphere

Materialize a BackgroundAtmosphere by converting O₃ functions to fields and converting constant gases to the grid's float type.

This is called internally by RadiativeTransferModel constructors.

materialize_surface_property(x, grid [, solar_position])

Convert a surface property (albedo, emissivity) to the form the radiative-transfer solver stores: a Number becomes a grid-eltype scalar and a Field passes through. Extend the three-argument form for property sources that must be resolved against the grid and the solar epoch (e.g. an observed-albedo dataset); it falls back to the two-argument form.

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_tendencies(
    microphysics,
    names::Tuple,
    ρ,
    ℳ,
    𝒰,
    constants
) -> Tuple

Compute the tendencies of names together, as a tuple in the order given.

The gridless counterpart of compute_microphysical_tendencies!. The default maps microphysical_tendency over names; schemes whose process rates are coupled across species override it to evaluate their bundle once and distribute it.

See also microphysical_tendency, compute_microphysical_tendencies!.

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.

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

Compute MoistureMassFractions from a microphysical state ℳ and scheme-dependent specific moisture $qᵛᵉ$.

The input $qᵛᵉ$ is the scheme-dependent specific moisture: vapor for non-equilibrium schemes, or equilibrium moisture ($qᵉ = qᵛ + qᶜˡ$) for saturation adjustment schemes.

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.

moisture_prognostic_name(_::Nothing) -> Symbol

Return the prognostic moisture field name as a Symbol for the given microphysics scheme.

The physical meaning of the prognostic moisture field depends on the scheme:

• Nothing / non-equilibrium: :ρqᵛ (true vapor density)
• SaturationAdjustment: :ρqᵉ (equilibrium moisture density, diagnostically partitioned)

moisture_specific_name(microphysics) -> Symbol

Return the specific (per-mass) moisture field name by stripping the ρ prefix from moisture_prognostic_name.

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.

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.

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

Return tuple() - Nothing microphysics has no prognostic variables.

set_to_mean!(ref::ExnerReferenceState, model)

Exner analogue of the ReferenceState method, for split-explicit CompressibleDynamics. Recompute the base exner_function/pressure/density by re-running the same discrete Exner column integration the constructor uses, with the horizontal-mean liquid-ice potential temperature and vapor mass fraction of the current model state. On height-coordinate grids the recomputed reference is horizontally uniform. Terrain-following model resets use their specialized constant-height mean and per-column integration path.

Unlike the anelastic ReferenceState method there is no rescale_densities option: the Exner reference is only the perturbation-form base state, not the prognostic density (ρᵈ), so changing it does not require rescaling the density-weighted prognostics.

set_to_mean!(reference_state, model; rescale_densities=false)

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

When rescale_densities=true, density-weighted prognostic fields (ρs, ρqᵗ, ρu, etc.) are rescaled by ρᵣ_new / ρᵣ_old so that the specific quantities (s, qᵗ, u, etc.) are unchanged. When false (default), the density-weighted fields are left as-is and only diagnostics are recomputed.

specific_humidity(model) -> Any

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

This always returns the actual vapor field $qᵛ$ from the microphysical fields, regardless of microphysics scheme.

specific_prognostic_moisture(model) -> Any

Return the prognostic specific moisture field for model.

This is $qᵛ$ for non-equilibrium schemes or $qᵉ$ for saturation adjustment schemes.

specific_prognostic_moisture_from_total(
    _::Nothing,
    qᵗ,
    ℳ
) -> Any

Convert total specific moisture $qᵗ$ to the scheme-dependent specific moisture $qᵛᵉ$ by subtracting the appropriate condensate from the microphysical state $ℳ$.

For non-equilibrium schemes, $qᵛᵉ = qᵛ = qᵗ - qˡ$ (subtract all condensate). For saturation adjustment schemes, $qᵛᵉ = qᵉ = qᵗ - qʳ$ (subtract only precipitation). For Nothing microphysics, $qᵛᵉ = qᵗ$ (all moisture is vapor).

This is used by parcel models that store total moisture $qᵗ$ as the prognostic variable, to produce the correct input for moisture_fractions.

specific_thermodynamic_field(formulation)

Return the specific (per unit mass) thermodynamic field the formulation evolves — what its advection operator reconstructs, as opposed to the density-weighted prognostic named by thermodynamic_density_name.

standard_ozone_profile(z) -> Any

An idealized climatological ozone volume mixing ratio (mol/mol) as a function of height z (m): a weak tropospheric background increasing toward the tropopause, blended into a Gaussian stratospheric layer peaking near 25 km. Keeps the stratospheric column near radiative balance in deep-column simulations — without ozone the upper column is far from radiative equilibrium and destabilizes when the spectral fluxes recompute. Not a substitute for an observed or model ozone climatology.

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 — the prognostic thermodynamic variable in coupling-density-weighted ("flux") form (ρθ, ρs, ρE).

The weighting density is the dynamics' coupling density (see dynamics_density): the reference density ρᵣ on the anelastic core and the prognostic dry-air density ρᵈ on the compressible core. The generic name (ρθ) is therefore ρᵈθ on CompressibleDynamics; the intensive variable is recovered as θ = ρθ / dynamics_density(dynamics).

thermodynamic_density_name(formulation)

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

total_density(
    i,
    j,
    k,
    dry_density,
    microphysics,
    moisture_density,
    microphysical_fields
) -> Any

Total air density $ρ = ρᵈ + ρᵗ$ at (i, j, k): the dry-air density dry_density plus the total_condensate_density $ρᵗ$. This is the diagnosed total mass density used where total mass enters the physics — the gravitational/buoyancy term and the equation of state.

total_density(dynamics)

Return the total air density ρ = ρᵈ + Σρˣ used by the thermodynamics, scalar advection, equation of state, and buoyancy. Defaults to dynamics_density — correct for formulations with a single density (e.g. the anelastic reference density). CompressibleDynamics overrides it with a diagnosed total-density field, distinct from the coupling density ρᵈ.

total_pressure(dynamics)

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

transport_velocities(
    model
) -> NamedTuple{(:u, :v, :w), <:Tuple{Any, Any, Any}}

Return the velocity tuple used for scalar advection transport.

For standard (non-terrain) models, this is model.velocities. For terrain-following coordinates, the vertical component is replaced by the contravariant vertical velocity $\tilde{w}$.

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;
    solver
) -> KernelFunctionOperation{_A, _B, _C, _D, T, K, D} where {_A, _B, _C, _D, T, K<:Breeze.AtmosphereModels.Diagnostics.DewpointTemperatureKernelFunction, D<:Tuple}

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 solver keyword argument (default SecantSolver(reltol=1e-4, abstol=0, maxiter=10)) controls the secant iteration; its convergence criterion compares the vapor pressure residual against the actual vapor pressure $pᵛ$.

Example

using Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid; microphysics=SaturationAdjustment())set!(model, θ=300, qᵗ=0.01)T⁺ = DewpointTemperature(model)# outputKernelFunctionOperation 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⁺)# output1×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.056, min=287.474, mean=288.266

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 Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300, qᵗ=0.01)θᵉ = EquivalentPotentialTemperature(model)Field(θᵉ)# output1×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.851, mean=326.006

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 Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300, qᵗ=0.01)θˡⁱ = LiquidIcePotentialTemperature(model)Field(θˡⁱ)# output1×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 Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300, qᵗ=0.01)θ = PotentialTemperature(model)Field(θ)# output1×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{_A, _B, _C, _D, T, K, D} where {_A, _B, _C, _D, T, K<:Breeze.AtmosphereModels.Diagnostics.SaturationSpecificHumidityKernelFunction, D<:Tuple}
SaturationSpecificHumidity(
    model,
    flavor_symbol
) -> KernelFunctionOperation{_A, _B, _C, _D, T, K, D} where {_A, _B, _C, _D, T, K<:Breeze.AtmosphereModels.Diagnostics.SaturationSpecificHumidityKernelFunction, D<:Tuple}

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's specific moisture field. This is equivalent to the :total_moisture flavor under saturated conditions with no condensate; or in other words, if the 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 Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300, qᵗ=0.01)θᵇ = StabilityEquivalentPotentialTemperature(model)Field(θᵇ)# output1×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.851, mean=326.006

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 $s$.

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

\[s = 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 $s$, or :density to return $ρ s$.

Examples

using Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300)s = StaticEnergy(model)Field(s)# output1×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.03019e5, min=302661.0, mean=3.0284e5

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 \frac{Rᵐ}{Rᵈ} = T \left( 1 + δᵛ qᵛ - qˡ - qⁱ \right)\]

The virtual potential temperature is defined analogously,

\[θᵛ = T \left( \frac{pˢᵗ}{p} \right)^{Rᵈ/cᵖᵈ} \left( 1 + δᵛ qᵛ - qˡ - qⁱ \right) .\]

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

using Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300, qᵗ=0.01)θᵛ = VirtualPotentialTemperature(model)Field(θᵛ)# output1×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.82, min=301.8, mean=301.81

azimuthal_mean!(profile, field; center, m) -> Any

Remap field (on an $(x, y, z)$ grid) onto the radial rings of profile (on an $(r, z)$ grid) about center, in place, by area-weighted binning of m × m sub-cells per Cartesian cell. profile and field must share their vertical grid; the radial rings are profile's uniform x-cells.

azimuthal_mean(field; radius, Nr, center, m)

Azimuthally average field into radial rings about center = (xc, yc), returning a Field on a one-dimensional-in-radius grid — Bounded in radius, Flat in azimuth, and Bounded in z — with Nr uniform rings spanning $[0, \texttt{radius}]$ and the same vertical grid as field.

Each Cartesian cell is split into an m × m block of sub-cells that are binned by radius, so a cell contributes to every ring it overlaps — uniform sub-sampling is area-weighting, which makes this a first-order conservative remap onto the radial rings (a pragmatic stand-in for a reduction on a true cylindrical grid). The kernel runs on the CPU and the GPU. Larger m resolves the rings more finely and keeps near-center rings populated; any ring that still catches nothing (only when $\texttt{radius}/N_r$ is finer than a sub-cell) is filled with NaN, not zero, so it reads as "no data" and doesn't bias a subsequent radial average.

using Oceananigans, Breezegrid = RectilinearGrid(size=(64, 64, 4), x=(-1, 1), y=(-1, 1), z=(0, 1),                       topology=(Periodic, Periodic, Bounded))c = CenterField(grid)set!(c, (x, y, z) -> 5)            # a constant fieldc̄ = azimuthal_mean(c; radius=1, Nr=8)maximum(c̄)                         # the azimuthal mean of a constant is that constant# output5.0

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

Create a bulk drag function for computing wall momentum fluxes using bulk aerodynamic formulas. The momentum flux is computed in the same form as the scalar bulk fluxes,

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

where Cᴰ is the drag coefficient, |U| = √(u² + v² + gustiness²) is the wind speed tangential to the wall (with gustiness to prevent singularities at low wind), u is the velocity component at the first cell face, and ρˢ is the surface density computed from the surface pressure and surface temperature. Monin–Obukhov similarity is a profile law for u (not ρu), so using u here keeps the formulation consistent with the similarity theory underlying Cᴰ.

The drag may be placed on any of the six boundaries of a bounded domain, on either of the two momentum components tangential to that wall: ρu and ρv on the bottom and top, ρv and ρw on the west and east, ρu and ρw on the south and north. The sign above is for the bottom; on every wall the drag removes tangential momentum from the domain.

When a FilteredSurfaceVelocities is supplied via filtered_velocities, every field entering the formula — the wind speed |U|, the velocity u, and the surface-layer virtual potential temperature difference Δθᵥ used in stability — is read from the filtered state. The surface density ρˢ is computed from the (slowly varying) surface temperature and pressure and is not filtered. Temporal filtering of the matching velocity is used to mitigate log-layer mismatch in wall-modeled large-eddy simulations, where the spurious correlation between the instantaneous friction velocity and matching-velocity fluctuations otherwise biases the surface stress (Nishizawa & Kitamura (2018); Shin, Yang & Howland (2025)). Filtering is supported on the bottom boundary only.

Monin–Obukhov consistency

ρˢ is computed from the surface temperature and the live model pressure extrapolated hydrostatically from the first cell center to the bottom face. It is therefore a true surface density that follows both terrain and changes in the model state.

Default surface temperature

If the user does not supply surface_temperature, materialization calls default_drag_surface_temperature(dynamics, …). The default exists for AnelasticDynamics (recovered from the reference state via Exner) but raises for CompressibleDynamics, which has no equivalent reference profile — pass surface_temperature explicitly in that case.

Keyword Arguments

• direction: The direction of the momentum component (XDirection(), YDirection(), or ZDirection()). 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, used to compute ρˢ and required when using PolynomialCoefficient with stability correction. Can be a Field, Function, or Number. A function takes the non-Flat coordinates of the wall followed by the time, as for Oceananigans boundary conditions: (x, y, t) on the bottom and top, (y, z, t) on the west and east, (x, z, t) on the south and north. (default: nothing)
• filtered_velocities: A FilteredSurfaceVelocities for temporally filtered wind speed, near-surface velocity, and θᵥ in the bulk formula. If nothing (default), instantaneous fields are used.

BulkSensibleHeatFluxFunction(
;
    coefficient,
    gustiness,
    surface_temperature,
    filtered_velocities
)

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 tangential to the wall, and $Δϕ$ is the difference between the near-wall atmospheric value and the wall value of the thermodynamic variable appropriate to the formulation:

• For LiquidIcePotentialTemperatureFormulation: $Δϕ = θ - θˢ$, where $θˢ = Tˢ / Πˢ$ and $Πˢ = (pˢ / pˢᵗ)^{Rᵈ / cᵖᵈ}$ (potential temperature flux)
• For StaticEnergyFormulation: $Δϕ = s - (cᵖᵐ Tˢ + g zˢ)$ (static energy flux), with $zˢ$ the height of the wall

Here $pˢ$ is the actual surface pressure, while $pˢᵗ$ is the fixed reference pressure used to define potential temperature.

The flux may be placed on any of the six boundaries of a bounded domain. The sign above is for the bottom; on every wall the flux carries heat into the domain when the wall is warmer than the adjacent air.

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 wall temperature. Can be a Field, a Function, or a Number. Functions are evaluated at the wall at every time step with the non-Flat coordinates of the wall followed by the time, as for Oceananigans boundary conditions: (x, y, t) on the bottom and top, (y, z, t) on the west and east, (x, z, t) on the south and north, and for example (x, t) on the bottom of a grid that is Flat in y.
• filtered_velocities: Either nothing (default) or FilteredSurfaceVelocities. Note that when filtered_velocities is not nothing, then automatically there is filtering in the scalar fields via FilteredSurfaceScalar with the same parameters (e.g., height, timescale) as filtered_velocities. Filtering is supported on the bottom boundary only.

BulkVaporFluxFunction(; coefficient, gustiness=0, surface_temperature,
                        surface_relative_humidity=1, moisture_availability=nothing,
                        filtered_velocities=nothing)

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

\[Jᵛ = - ρˢ Cᵛ |U| (qᵛ - qˢ), \qquad qˢ = β ℋˢ qᵛ⁺(Tˢ) + (1 - β) qᵛ,\]

where $Cᵛ$ is the transfer coefficient, $|U|$ is the wind speed tangential to the wall, $qᵛ$ is the near-wall specific humidity, and $qˢ$ is the specific humidity of the air in contact with the wall. Over the wet fraction $β$ of the wall (the moisture_availability) that is the saturation specific humidity $qᵛ⁺$ at the wall temperature $Tˢ$ times the wall relative humidity $ℋˢ$ (unity for a wet wall); over the dry fraction it is the humidity of the air itself, so that $qᵛ - qˢ = β (qᵛ - ℋˢ qᵛ⁺)$ and the flux is $β$ times the flux over a wet wall.

The flux may be placed on any of the six boundaries of a bounded domain. The sign above is for the bottom; on every wall the flux carries vapor into the domain when the wall is moister than the adjacent air.

Keyword Arguments

• coefficient: The vapor transfer coefficient.
• gustiness: Minimum wind speed to prevent singularities (default: 0).
• surface_temperature: The wall temperature. Can be a Field, a Function, or a Number. Used to compute the saturation specific humidity at the wall. Functions take the non-Flat coordinates of the wall followed by the time, as for Oceananigans boundary conditions.
• surface_relative_humidity: The relative humidity of the air in contact with the wall, between 0 and 1 (default: 1, a saturated wall). Can be a Field, a Function, or a Number.
• moisture_availability: The fraction $β ∈ [0, 1]$ of the wall that is wet. nothing (default) takes the value carried by a PolynomialCoefficient coefficient, whose stability correction uses the same surface humidity, and 1 (a wet wall, an ocean) for a constant coefficient. A value that disagrees with a PolynomialCoefficient is an error. The phase of the surface water follows the coefficient in the same way, and is liquid for a constant coefficient.
• filtered_velocities: Either nothing (default) or FilteredSurfaceVelocities. Note that when filtered_velocities is not nothing, then automatically there is filtering in the scalar fields via FilteredSurfaceScalar with the same parameters (e.g., height, timescale) as filtered_velocities. Filtering is supported on the bottom boundary only.

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ᵖᵐ$ and the Exner function $Π$.

The relationship is:

\[Jᶿ = 𝒬ᵀ / (cᵖᵐ Π)\]

where $𝒬ᵀ$ is the energy flux and $Jᶿ$ is the potential temperature flux. At fixed moisture $T = Π θ + (ℒˡᵣ qˡ + ℒⁱᵣ qⁱ) / cᵖᵐ$ gives $δT = Π δθ$, so an enthalpy input $𝒬ᵀ = ρ cᵖᵐ \overline{w'T'}$ reaches $θ$ as $Jᶿ = ρ \overline{w'T'} / Π = 𝒬ᵀ / (cᵖᵐ Π)$. The condensate term cancels, so the conversion holds saturated as well as dry. This is the same conversion the $ρE$ interior forcing and the radiative flux divergence apply. In the same spirit, BulkSensibleHeatFlux refers its wall temperature to $θ$, though through a dry $Πˢ$ at the wall face rather than the moist $Π$ of the adjacent cell.

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

FilteredSurfaceScalar(grid; height=nothing, filter_timescale=Inf)

A two-dimensional field storing a temporally filtered near-surface scalar for use in bulk flux boundary conditions.

The filter update is the same exponential form as FilteredSurfaceVelocities.

Keyword Arguments

• height: Reference height (m) for scalar evaluation. If nothing, the first grid cell center value is used.
• filter_timescale: Filter time scale τ in seconds (default: Inf).

FilteredSurfaceVelocities(grid; height=nothing, filter_timescale=Inf)

Two-dimensional fields storing temporally filtered near-surface velocities and the surface-layer virtual potential temperature difference for use in bulk flux boundary conditions. Filtering the matching velocity mitigates log-layer mismatch in wall-modeled large-eddy simulations by removing the spurious correlation between the instantaneous friction velocity and matching-velocity fluctuations (Nishizawa & Kitamura (2018); Shin, Yang & Howland (2025)).

The filtered velocities ū, v̄ (and the surface-layer difference Δθ̄ᵥ when a stability-dependent bulk coefficient is attached) are updated each time step via an exponential (first-order) filter:

ū ← (ū + ϵ u_new) / (1 + ϵ),     ϵ = Δt / τ

where τ is the filter_timescale.

Δθ̄ᵥ filters the result the stability correction needs, $Δθᵥ = θᵥ(z₁) - θᵥˢ$, formed at every update by the attached PolynomialCoefficient from the instantaneous first-cell virtual potential temperature and specific humidity, the surface temperature, the moisture availability and the surface phase. Filtering the one difference rather than each of its inputs stores a single field, and is exact since the filter is linear. When the bulk coefficient carries no stability correction the field is allocated but unused.

Keyword Arguments

• height: Reference height (m) for velocity evaluation. If nothing (default), the first grid cell center value is used. If a number, velocity is linearly interpolated to that height. (Δθ̄ᵥ is always formed at the first cell center, matching the height at which bulk_coefficient evaluates stability.)
• filter_timescale: Filter time scale τ in seconds (default: Inf, no filtering).

FittedStabilityFunction(scalar_roughness_length;
                        richardson_number_mapping = RichardsonNumberMapping(typeof(scalar_roughness_length)),
                        stability_function = StabilityFunction(typeof(scalar_roughness_length)))

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: StabilityFunction (default: Hogström (1996) / Beljaars & Holtslag (1991)).

References

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

PolynomialCoefficient(
;
    ...
) -> PolynomialCoefficient{_A, Nothing, SF, PlanarLiquidSurface, Nothing, Nothing, Nothing, Nothing} where {_A, SF<:(FittedStabilityFunction{_A, RM, SP} where {_A, RM<:Breeze.BoundaryConditions.RichardsonNumberMapping, SP<:Breeze.BoundaryConditions.StabilityFunction})}
PolynomialCoefficient(
    FT;
    polynomial,
    roughness_length,
    minimum_wind_speed,
    stability_function,
    surface,
    moisture_availability,
    transfer_type
) -> PolynomialCoefficient{_A, Nothing, SF, PlanarLiquidSurface, Nothing, Nothing, Nothing, Nothing} where {_A, SF<:(FittedStabilityFunction{_A, RM, SP} where {_A, RM<:Breeze.BoundaryConditions.RichardsonNumberMapping, SP<:Breeze.BoundaryConditions.StabilityFunction})}

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 and Yeager (2009) form:

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

where $U_h$ is the wind speed at measurement height $h$. The polynomial evaluates to the coefficient itself; Large and Yeager publish their coefficients in units of $10^{-3}$, a scaling the default polynomials carry explicitly.

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 and Yeager (2009) polynomial will be automatically selected based on the boundary condition type:

• BulkDrag: default_neutral_drag_polynomial = (0.142, 0.076, 2.7) .* 1e-3 for momentum
• BulkSensibleHeatFlux: default_neutral_sensible_heat_polynomial = (0.128, 0.068, 2.43) .* 1e-3 for sensible heat
• BulkVaporFlux: default_neutral_latent_heat_polynomial = (0.120, 0.070, 2.55) .* 1e-3 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: The phase of the surface water, which selects the saturation specific humidity at the surface in the stability correction: PlanarLiquidSurface() (default), PlanarIceSurface() or PlanarMixedPhaseSurface(liquid_fraction).
• moisture_availability: The fraction $β ∈ [0, 1]$ of the surface that is saturated (default: 1, an ocean). The surface specific humidity entering the stability correction is $qˢ = β qᵛ⁺(Tˢ) + (1 - β) qᵛ$, with $qᵛ$ the specific humidity of the air in the first cell, so that $β = 0$ describes a dry surface whose virtual potential temperature carries no moisture contribution of its own. See wall_virtual_potential_temperature.

The measurement height is automatically determined from the grid as half the first-cell thickness, the height of its center above the local surface.

Examples

using Breeze.BoundaryConditions: PolynomialCoefficient# Polynomial coefficient with default settingscoef = PolynomialCoefficient()# outputPolynomialCoefficient{Float64}├── polynomial: nothing├── roughness_length: 0.00015 m├── minimum_wind_speed: 0.1 m/s├── surface: PlanarLiquidSurface├── moisture_availability: 1.0└── stability_function: FittedStabilityFunction (Li et al. 2010)

using Breeze.BoundaryConditions: PolynomialCoefficient# With explicit polynomialcoef = PolynomialCoefficient(polynomial = (0.000142, 7.6e-5, 0.0027))# outputPolynomialCoefficient{Float64}├── polynomial: (0.000142, 7.6e-5, 0.0027)├── roughness_length: 0.00015 m├── minimum_wind_speed: 0.1 m/s├── surface: PlanarLiquidSurface├── moisture_availability: 1.0└── stability_function: FittedStabilityFunction (Li et al. 2010)

using Breeze.BoundaryConditions: PolynomialCoefficient# No stability correctioncoef = PolynomialCoefficient(stability_function = nothing)# outputPolynomialCoefficient{Float64}├── polynomial: nothing├── roughness_length: 0.00015 m├── minimum_wind_speed: 0.1 m/s├── surface: PlanarLiquidSurface├── moisture_availability: 1.0└── stability_function: Nothing

References

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

RichardsonNumberMapping(FT = Oceananigans.defaults.FloatType;
                        stable_unstable_transition = 0,
                        strongly_stable_transition = 0.2,
                        aᵘ₁₁ =  0.0450, bᵘ₁₁ =  0.0030, bᵘ₁₂ =  0.0059,
                        aᵘ₂₁ = -0.0828, aᵘ₂₂ =  0.8845,
                        bᵘ₃₁ =  0.1739, bᵘ₃₂ = -0.9213, bᵘ₃₃ = -0.1057,
                        aʷ₁₁ =  0.5738, aʷ₁₂ = -0.4399,
                        aʷ₂₁ = -4.901,  aʷ₂₂ = 52.50,
                        bʷ₁₁ = -0.0539, bʷ₁₂ =  1.540,
                        bʷ₂₁ = -0.6690, bʷ₂₂ = -3.282,
                        aˢ₁₁ =  0.7529, aˢ₂₁ = 14.94,
                        bˢ₁₁ =  0.1569, bˢ₂₁ = -0.3091, bˢ₂₂ = -1.303)

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.

StabilityFunction(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

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

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ᵖᵐ$ and the Exner function $Π$.

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 wall momentum drag, on any of the six boundaries.

See BulkDragFunction for details.

Examples

using Breezedrag = BulkDrag(coefficient=1e-3, gustiness=0.1)# outputFluxBoundaryCondition: BulkDragFunction(direction=Nothing, coefficient=0.001, gustiness=0.1)

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

using Oceananigans.Grids: XDirectionu_drag = BulkDrag(direction=XDirection(), coefficient=1e-3)ρu_bcs = FieldBoundaryConditions(bottom=u_drag)# outputOceananigans.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. The same condition may be placed on the walls of a closed box; the direction is inferred from the momentum component it is attached to:

drag = BulkDrag(coefficient=1e-3)ρv_bcs = FieldBoundaryConditions(west=drag, east=drag, bottom=drag, top=drag)# outputOceananigans.FieldBoundaryConditions, with boundary conditions├── west: FluxBoundaryCondition: BulkDragFunction(direction=Nothing, coefficient=0.001, gustiness=0)├── east: FluxBoundaryCondition: BulkDragFunction(direction=Nothing, coefficient=0.001, gustiness=0)├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)├── bottom: FluxBoundaryCondition: BulkDragFunction(direction=Nothing, coefficient=0.001, gustiness=0)├── top: FluxBoundaryCondition: BulkDragFunction(direction=Nothing, coefficient=0.001, gustiness=0)└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

BulkSensibleHeatFlux(; coefficient, gustiness=0, surface_temperature)

Create a FluxBoundaryCondition for wall sensible heat flux, on any of the six boundaries.

The bulk formula computes

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

where $Δϕ$ depends on the thermodynamic formulation: $Δθ$ for potential temperature or $Δs$ for static energy. The formulation is set automatically during model construction.

See BulkSensibleHeatFluxFunction for details.

Example

using BreezeTˢ(x, y, t) = 290 + 2 * sign(cos(2π * x / 20e3))ρE_bc = BulkSensibleHeatFlux(coefficient = 1e-3,                             gustiness = 0.1,                             surface_temperature = Tˢ)# outputFluxBoundaryCondition: BulkSensibleHeatFluxFunction(coefficient=0.001, gustiness=0.1)

BulkVaporFlux(; coefficient, surface_temperature, surface_relative_humidity=1,
                moisture_availability=nothing, gustiness=0)

Create a FluxBoundaryCondition for wall moisture flux, on any of the six boundaries.

The specific humidity of the air in contact with the wall is computed from surface_temperature and surface_relative_humidity (unity by default, a wet wall); moisture_availability is the fraction of the wall that is wet, 1 by default.

See BulkVaporFluxFunction for details.

Example

using BreezeTˢ(x, y, t) = 290 + 2 * sign(cos(2π * x / 20e3))moisture_bc = BulkVaporFlux(coefficient = 1e-3,                            gustiness = 0.1,                            surface_temperature = Tˢ)# outputFluxBoundaryCondition: 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ᵖᵐ$ and the Exner function $Π$ 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ᵖᵐ$ and the Exner function $Π$ to obtain the energy flux: $𝒬ᵀ = Jᶿ cᵖᵐ Π$.

cos_solar_zenith_angle(
    i,
    j,
    grid::RectilinearGrid{<:Any, <:Flat, <:Flat, <:Bounded},
    datetime::Dates.AbstractDateTime
) -> 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.AbstractDateTime,
    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.AbstractDateTime) -> Int64

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

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.AbstractDateTime,
    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.

abstract type AbstractRamp

Abstract supertype for upper-sponge ramp shapes. A concrete AbstractRamp is callable as (ramp)(z, sponge_top, depth) and returns a value in $[0, 1]$: zero below $z_{\rm sponge\_top} - \text{depth}$, rising to one at the lid $z = z_{\rm sponge\_top}$.

abstract type AcousticDampingStrategy

Abstract supertype for divergence damping applied inside the substep loop.

Concrete subtypes:

• NoDivergenceDamping — no damping.
• ThermalDivergenceDamping — Klemp, Skamarock & Ha (2018) momentum correction using the discrete $δ_τ(ρθ)$ tendency as the divergence proxy. This is the default used by SplitExplicitTimeDiscretization.

abstract type AcousticOuterScheme

Abstract supertype for the outer Runge–Kutta scheme that drives the acoustic substep loop. The current implementation supports a single concrete subtype:

• WickerSkamarock3 — three-stage Wicker–Skamarock RK3 with stage fractions $β = (1/3, 1/2, 1)$. This is the only outer scheme supported today and the default for SplitExplicitTimeDiscretization.

The interface exists to make the outer-scheme commitment explicit in the type system and to provide a clean extension point for a future Multirate Infinitesimal Step (MIS) outer scheme. A concrete subtype is expected to provide a stage_fractions method returning its stage-fraction tuple.

abstract type AcousticSubstepDistribution

Abstract supertype for the choice of how acoustic substeps are distributed across the three Wicker–Skamarock RK3 stages.

Concrete subtypes:

• ProportionalSubsteps — each stage independently covers its own interval $β Δt$ with $Nτ = ⌈β N⌉$ substeps of size $Δτ = β Δt / Nτ$ (count proportional to the stage fraction; size fitted so the substeps exactly tile $β Δt$). This is the default.

• ConstantSubstepSize — every stage uses the same substep size $Δτ = Δt/N$ ($N$ rounded up to a multiple of 6 so $β N$ is integral), with stage-dependent counts $Nτ = β N$.

• MonolithicFirstStage — stage 1 collapses to a single substep of size $Δt/3$; stages 2 and 3 are the same as ConstantSubstepSize.

struct AcousticSubstepper{N, FT, D, AD, US, CF, MP, TAV, GT, TS, WC, DC, TWC}

Storage and parameters for the split-explicit acoustic substepper (scheme described in the module header). Πᴸ=(pᴸ/pˢᵗ)^κ, θᴸ=ρθᴸ/ρᴸ, γᵐRᵐᴸ are cached once per stage (recomputing inline per call is much slower on H100); ρᴸ, ρθᴸ, pᴸ and the stage-entry momenta are read live from model.dynamics.* / model.momentum.* (untouched by the loop) and are the recovery base for _recover_full_state! — no snapshot fields. The vertical solve is a (possibly off-centered) Crank-Nicolson tridiagonal Schur system for (ρw)′.

Fields:

• substeps: acoustic substeps N per Δt (nothing ⇒ adaptive via acoustic_cfl).
• acoustic_cfl: target horizontal acoustic Courant number for the adaptive count (default 0.5).
• forward_weight: CN off-centering ω (0.5 = centered; default 0.65).
• damping, substep_distribution: divergence-damping strategy; substep allocation across WS-RK3 stages.
• linearization_exner (Πᴸ), linearization_potential_temperature (θᴸ), linearization_gamma_R_mixture (γᵐRᵐᴸ, the moist PGF coefficient): per-stage caches.
• density_perturbation (ρ′), density_potential_temperature_perturbation ((ρθ)′), momentum_perturbation ((ρu/v/w)′ as .u/.v/.w): perturbation prognostics advanced in the loop.
• density_predictor, density_potential_temperature_predictor: explicit predictors before the vertical solve.
• previous_density_potential_temperature_perturbation: prior-substep (ρθ)′, for Klemp 2018 damping.
• time_averaged_velocities: acoustic-mean velocity for non-acoustic scalar transport (moisture/tracers/ chemistry/TKE); the slow ρθ tendency uses the current RK predictor velocity instead, not this cache.
• slow_vertical_momentum_tendency (Gˢρw, z-faces): advection+Coriolis+closure+forcing (PGF/buoyancy excluded — those are in the fast operator).
• vertical_solver_source_term (z-faces): explicit RHS of the (ρw)′ tridiagonal system.
• vertical_solver: BatchedTridiagonalSolver for the implicit (ρw)′ update.
• vertical_velocity_cache, density_cache: the stage-entry predictor velocity and its carrier dry density, frozen for implicit_substep! (see cache_advecting_state!); nothing without adaptive-implicit advection.
• time_averaged_vertical_velocity_cache: the acoustic-mean transport velocity the moisture and tracer tendencies were built with, frozen for scalar_substep! (see cache_transport_velocity!); nothing without adaptive-implicit advection.

AcousticSubstepper(
    grid,
    split_explicit::SplitExplicitTimeDiscretization;
    prognostic_momentum,
    substep_floattype,
    cache_advecting_state
) -> AcousticSubstepper{_A, _B, _C, AD, _D, CF, MP, TAV, GT, TS, WC, DC, TWC} where {_A, _B, _C, AD<:AcousticSubstepDistribution, _D, CF<:(Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), MP<:(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}}}), TAV<:(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}}}), GT<:(Field{Center, Center, Face, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}), TS<:(Oceananigans.Solvers.BatchedTridiagonalSolver{Breeze.CompressibleEquations.AcousticTridiagLower, Breeze.CompressibleEquations.AcousticTridiagDiagonal, Breeze.CompressibleEquations.AcousticTridiagUpper, _A, G, Nothing, Oceananigans.Grids.ZDirection} where {_A, G<:Oceananigans.Grids.AbstractGrid}), WC<:Union{Nothing, Field{Center, Center, Face, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}, DC<:Union{Nothing, Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}, TWC<:Union{Nothing, Field{Center, Center, Face, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}}

Construct an AcousticSubstepper. The perturbation face fields $(ρu)′, (ρv)′, (ρw)′$ and the scalar-transport velocities use topology-derived BCs (periodic wrap / impenetrability), not the prognostic momentum's BCs: inheriting them would imprint the full-state wall target onto the perturbation halo for a nonzero NormalFlowBoundaryCondition (issue #716) and apply momentum BCs to velocity fields. The wall target re-enters via the prognostic momentum's own BC after each substep's momentum update. The prognostic_momentum kwarg is retained for backwards compatibility but no longer consulted.

struct CompressibleDynamics{TD, D, DT, P, FT, RS, TM, CV, CM}

Fully compressible dynamics with prognostic density and diagnostic pressure.

Fields

• dry_density: Prognostic dry-air density field ρᵈ
• total_density: Diagnosed total air density ρ = ρᵈ + Σρˣ (used for thermodynamics, scalar advection, EOS, buoyancy)
• pressure: Diagnostic pressure field p = ρ Rᵐ T
• standard_pressure: Reference pressure pˢᵗ for potential temperature (default 10⁵ Pa)
• base_pressure: Mean pressure of the reference atmosphere at $z = 0$ (p₀). The pressure at the ground is surface_pressure, which the reference state derives from it
• time_discretization: Time discretization scheme (SplitExplicitTimeDiscretization or ExplicitTimeStepping)
• reference_state: The single fixed hydrostatically-balanced reference state for base-state pressure/buoyancy correction (perturbation-form PGF), or nothing when disabled. An ExnerReferenceState whose fields are grid-polymorphic: a 1D column on height-coordinate grids, and horizontally-varying 3D fields on terrain-following grids (where a single column is not hydrostatically consistent per terrain column).
• terrain_metrics: TerrainMetrics for terrain-following coordinates (or nothing). This — not reference_state — is the sole "is this a terrain grid?" signal.
• w̃, ρw̃: contravariant vertical velocity / momentum diagnostic fields (or nothing when no terrain metrics)

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)

The moist equation-of-state θˡⁱ→T temperature inversion is controlled by the thermodynamic formulation, not the dynamics: see temperature_solver on LiquidIcePotentialTemperatureFormulation.

CompressibleDynamics(
;
    ...
) -> CompressibleDynamics{ExplicitTimeStepping, Nothing, Nothing, Nothing, Float64, Breeze.CompressibleEquations.AutoReference, SlopeOutsideInterpolation, Nothing, Nothing}
CompressibleDynamics(
    time_discretization;
    standard_pressure,
    base_pressure,
    reference_potential_temperature,
    reference_temperature,
    reference_vapor_mass_fraction,
    slope_stencil,
    terrain_metrics,
    reference_state,
    temperature_tolerance,
    temperature_maxiter,
    surface_pressure
) -> CompressibleDynamics{_A, Nothing, Nothing, Nothing, Float64, Breeze.CompressibleEquations.AutoReference, SlopeOutsideInterpolation, Nothing, 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)

• base_pressure: Mean pressure of the reference atmosphere at $z = 0$ (default: 101325.0 Pa). A datum, not the pressure at the ground: over terrain, or on a domain whose bottom is raised, the two differ by $O(ρgh)$

• 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, which uses the automatic θᵣ = 288 K profile when reference_state = :auto. When provided, this profile replaces the automatic profile when building the ExnerReferenceState.

• reference_vapor_mass_fraction: Optional vapor mass fraction for building a moist compressible reference state. Can be a constant qᵛ, function qᵛ(z), or field, and is used with reference_potential_temperature.

• slope_stencil: Pressure-gradient slope-interpolation stencil for terrain-following grids. Default: SlopeOutsideInterpolation. Ignored on non-terrain-following grids.

• terrain_metrics: Escape hatch — pass a pre-built TerrainMetrics to bypass the automatic build. Default: nothing (auto-build from the grid using slope_stencil).

• reference_state: Whether to carry the single hydrostatic reference state used for the perturbation-form pressure-gradient force and buoyancy. Default: :auto — on a bounded vertical grid, build a standard-atmosphere (θᵣ = 288 K) hydrostatic reference: a 1D column on height-coordinate grids, 3D fields on terrain-following grids. Periodic and flat vertical topologies carry no automatic reference because a nontrivial hydrostatic atmosphere is incompatible with periodicity and unnecessary without a vertical dimension. Pass reference_state = nothing to disable it entirely — the PGF and buoyancy then difference the full pressure, reproducing the un-corrected behavior (useful for testing). Disabling is mutually exclusive with an explicit reference profile. To replace the reference with one deduced from an initial state's horizontal mean, call set!(model; …, compute_reference_state=true).

  Deep near-isentropic reference profiles

  The reference integrates the hydrostatic equation up each column using θᵣ(z). A (nearly) constant-θ column is isentropic and its hydrostatic pressure reaches zero at a finite height (≈ cᵖ θ / g, about 29 km for the default θᵣ = 288 K); if the domain top exceeds that height the integration has no positive-pressure solution and the reference fills with NaN. Physical, stably-stratified profiles are unaffected. For such a deep, near-isentropic setup pass either a stratified reference_potential_temperature or reference_state = nothing (full-pressure form).

struct ConstantSubstepSize <: AcousticSubstepDistribution

Acoustic substep distribution where every stage uses the same substep size $Δτ = Δt/N$. $N$ is rounded up to a multiple of 6 (= LCM of the WS-RK3 stage denominators 2 and 3) so the per-stage count $Nτ = β_\mathrm{stage} N$ is an exact integer and each stage covers exactly $β Δt$ — uniform Δτ, at the cost of over-resolving (substep count is the next multiple of 6 ≥ the CFL minimum).

struct CubicRamp <: AbstractRamp

Hermite cubic "smoothstep" sponge ramp. $s² (3 − 2s)$ where $s = \text{clamp}((z − (H − \text{depth}))/\text{depth}, 0, 1)$.

Has zero derivative at both the layer base and the lid, so absorbs upgoing waves smoothly without the reflective kink of LinearRamp. Functionally equivalent to Sin2Ramp but ~5–10× cheaper inside the GPU kernel (no transcendental). Recommended default.

struct DirectDivergenceDamping{FT} <: AcousticDampingStrategy

Acoustic divergence damping that forms the horizontal θ-flux divergence $δ = ∂ₓ(θᴸ(ρu)′) + ∂_y(θᴸ(ρv)′)$ directly from the perturbation momentum, rather than approximating it through the $(ρθ)′$ substep tendency the way ThermalDivergenceDamping does (Klemp, Skamarock & Ha 2018, their eq. 36). After each acoustic substep the horizontal perturbation momentum receives the correction

\[Δ(ρu)′ = α\, Δx²\, ∂ₓ δ / θᴸ, \qquad Δ(ρv)′ = α\, Δy²\, ∂_y δ / θᴸ,\]

with the single dimensionless coefficient α (MPAS config_smdiv, default 0.1; the Laplacian-diffusion stability bound is α ≲ 0.2). The divergence is horizontal only: the damped quantity must match the divergence in the $Θ = ρθ$ equation, and folding in the vertical θ-flux divergence damps the resolved vertical flux and destabilizes the flow. Differencing the velocity field directly (rather than the $(ρθ)′$ tendency) carries no $1/Δτ$ in the diffusivity, which also avoids the thermal proxy's cold-start $∝ α/Δτ$ spurious force (cf. PR #794).

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 LinearRamp <: AbstractRamp

Linear sponge ramp. Cheap but introduces a kink at the bottom of the sponge layer (nonzero slope at $z = H − \text{depth}$), which can cause small partial reflection of upgoing waves in idealised tests. WRF's older damp_opt = 2 form uses this.

struct MonolithicFirstStage <: AcousticSubstepDistribution

Acoustic substep distribution where stage 1 collapses to a single substep of size $Δt/3$; stages 2 and 3 are the same as ConstantSubstepSize ($N/2$ and $N$ substeps of size $Δτ = Δt/N$).

struct NoDivergenceDamping <: AcousticDampingStrategy

No acoustic divergence damping. The substep loop advances perturbation fields without applying any post-substep momentum correction.

struct ProportionalSubsteps <: AcousticSubstepDistribution

Acoustic substep distribution where each WS-RK3 stage independently covers its interval $β_\mathrm{stage} Δt$ with $Nτ = ⌈β_\mathrm{stage} N⌉$ substeps of size $Δτ = β_\mathrm{stage} Δt / Nτ$. The count is proportional to the stage fraction and the size is fitted so the substeps exactly tile each stage — exact coverage at the minimum substep count (no global quantization; Δτ may differ slightly by stage).

This is the default.

struct Sin2Ramp <: AbstractRamp

$\sin^2$ sponge ramp from Klemp, Dudhia & Hassiotis (2008). Same zero-derivative-at-both-ends behaviour as CubicRamp, but with a transcendental call. Provided for parity with WRF (damp_opt = 3) / MPAS-Atmosphere; prefer CubicRamp for performance in new code.

struct SplitExplicitTimeDiscretization{N, FT, D, US, AD<:AcousticSubstepDistribution}

Time discretization for fully compressible dynamics that integrates slow terms with a Wicker-Skamarock RK3 outer loop and acoustic terms with split-explicit inner substeps.

The constructor accepts substeps or an acoustic_cfl for choosing the number of acoustic substeps, a forward_weight for off-centering the acoustic solve, an acoustic damping strategy such as ThermalDivergenceDamping, an optional UpperSponge, and a substep_distribution such as ProportionalSubsteps.

Backward integration (time_step!(model, Δt) with Δt < 0) is supported for the linearized acoustic substep loop. See the field-documentation docstring for the A-stability argument, sign-handling of the adaptive substep count, and the irreversibility caveat for the optional UpperSponge.

struct ThermalDivergenceDamping{FT, LS} <: AcousticDampingStrategy

Acoustic divergence damping that uses the (ρθ)′ tendency as a discrete proxy for the momentum divergence. From the linearized ρθ-continuity equation $\partial_t (ρθ)' + \nabla\cdot(ρθ^L u') = 0$, the per-substep quantity

\[D \equiv \frac{(ρθ)' - (ρθ)'_\mathrm{old}}{θ^L} \approx -Δτ \, \nabla\cdot(ρu)'\]

is what would otherwise require an extra divergence operator and an extra kernel pass. Building the correction from D reuses the substep's already-resident (ρθ)′ snapshots — that's the algorithmic choice this damping is named for.

Used by Klemp, Skamarock & Ha (2018) / Skamarock & Klemp (1992) / Baldauf (2010). After each acoustic substep, the horizontal momentum perturbation components $(ρu)′$ and $(ρv)′$ pick up an explicit correction proportional to the horizontal gradient of $D$. If damp_vertical = true, the vertical component is folded implicitly into the column tridiag as a Laplacian on the acoustic vertical momentum perturbation: $(ρw)′$ for height-coordinate dynamics and $(ρ ilde{w})′$ for terrain-following dynamics. By default, damp_vertical = false and vertical acoustic damping comes from the off-centered implicit solve.

Per-substep momentum correction (Klemp, Skamarock & Ha 2018 eq. 36, MPAS form):

\[Δ(ρu)′ = -γ · ∂_x D , \quad Δ(ρv)′ = -γ · ∂_y D .\]

with local per-direction horizontal diffusivities. On a uniform square grid this is the finite-difference analogue of MPAS's coef_divdamp = 2·smdiv·config_len_disp/Δτ:

\[γ_x = α \, Δx^2 / Δτ , \qquad γ_y = α \, Δy^2 / Δτ .\]

On anisotropic or latitude-longitude grids, using the local per-direction spacing keeps the nondimensional explicit damping strength approximately uniform across the mesh. Pass length_scale = ℓ to override the automatic local scale with the fixed diffusivity $γ = α ℓ^2 / Δτ$ when a nominal mesh length is more appropriate. The optional vertical tridiag contribution uses $γ_z = α Δz² / Δτ$ when damp_vertical = true.

$α$ is the dimensionless Klemp 2018 coefficient (= MPAS config_smdiv, default 0.1). The combined 2-D horizontal explicit-time stability bound is $8α ≤ 2 → α ≤ 0.25$; the default sits well below it. Combined with the SplitExplicitTimeDiscretization default $\omega = 0.65$, this preserves the exact discrete rest atmosphere at Δt = 20 s and damps divergent acoustic noise in production runs. It should not be read as a guarantee that every grid-scale balanced-mode growth diagnostic is physically correct.

Fields

• coefficient: Dimensionless damping coefficient $α$ (Klemp 2018 / MPAS config_smdiv). Default 0.1. The horizontal part is explicit and obeys the usual 8α ≤ 2 2-D combined CFL. When damp_vertical = true, the vertical contribution is implicit and is folded into the column tridiag.
• length_scale: Optional override for the dispersion length $d$. Default nothing (auto: local $γ_x = α Δx^2 / Δτ$ and $γ_y = α Δy^2 / Δτ$). Setting length_scale = ℓ forces a fixed $γ = α \, ℓ² / Δτ$ in both horizontal directions.
• damp_vertical: If true, the vertical part of the divergence damping is folded into the column tridiag (a Laplacian on (ρw)′ in height coordinates or (ρw̃)′ in terrain-following coordinates). If false (default), no extra vertical damping is applied — the vertical acoustic modes are damped solely by the off-centering of the implicit pressure-gradient solve ($\omega > 0.5$), which Klemp et al. 2018 eq. (32) shows is algebraically equivalent to a vertical divergence damping with diffusivity $γ_z = c² Δτ s/2$ where $s = 2\omega - 1$.

struct UpperSponge{FT, R<:AbstractRamp}

Implicit upper Rayleigh sponge for the substepper inner loop. Damps the acoustic vertical momentum perturbation toward zero inside a layer of thickness depth below the model lid, with peak damping rate damping_rate (in 1/s) at the lid scaled by ramp(z). The damped variable is $(ρw)′$ for height-coordinate dynamics and $(ρ ilde{w})′$ for terrain-following dynamics.

The damping is applied inside the column tridiag as a CN-weighted contribution (paralleling the existing implicit divergence-damping treatment): $δτᵐ⁺ × \text{rate} × \text{ramp}(z)$ on the LHS diagonal, $δτˢ⁻ × \text{rate} × \text{ramp}(z)$ on the explicit-half RHS. This matches the Rayleigh-layer form of the Klemp, Dudhia & Hassiotis (2008) absorbing treatment used in WRF (damp_opt=3) and MPAS-Atmosphere. The profile shape is controlled by ramp; use Sin2Ramp for the classic $\sin^2$ profile.

Keyword arguments

• damping_rate: peak damping rate at the lid, in 1/s. Default 0.2.

  Because the damping is fully implicit in the inner-loop tridiag, it is unconditionally stable for any positive value, so the choice is guided by physics rather than CFL. Typical guidance:

  • $\text{rate} ≳ N$ (Brunt–Väisälä frequency, ~0.01 /s in the stratosphere) is the lower bound at which gravity waves are absorbed rather than reflected within the layer crossing time.
  • WRF's dampcoef default and Klemp et al.'s recommendation is 0.2 (i.e. τ ≈ 5 s at the lid) — comfortably above $N$ and aggressive enough to absorb in one or two crossings.
  • Larger values are fine numerically but produce a sharper "cap" near the lid; if the application cares about resolved dynamics just below the sponge, prefer $\text{rate} ≈ 0.1$ and a deeper layer.

• depth: sponge-layer thickness below the lid, in metres along the reference vertical coordinate. Default 5e3.

  Should span at least ~10 grid cells in the vertical to give the smooth profile room to absorb without aliasing; for $Δz ≈ 1\,\text{km}$ the default of 5 km gives 5 cells (marginal — bump to 10 km if w-spectrum has structure right below the lid).

• ramp: an AbstractRamp controlling the profile shape. Default CubicRamp(). Other built-ins: Sin2Ramp(), LinearRamp(). Custom shapes are supported by subtyping AbstractRamp and defining (::MyRamp)(z, sponge_top, depth).

The ramp depends only on the reference vertical coordinate (no horizontal variation), so the sponge does not break zonal symmetry and remains uniform over terrain-following grids.

struct WickerSkamarock3 <: AcousticOuterScheme

Three-stage Wicker–Skamarock RK3 outer scheme (Wicker and Skamarock 2002) with canonical stage fractions $β = (1/3, 1/2, 1)$. Each stage resets the prognostic state to $U^n$ and applies a fraction $β_k Δt$ of the slow tendency evaluated at the previous-stage state, while the acoustic substep loop advances linearized perturbations about each RK stage-entry state.

acoustic_rk3_substep_loop!(
    model::AtmosphereModel,
    substepper,
    Δt,
    β_stage,
    Uᴸ
)
acoustic_rk3_substep_loop!(
    model::AtmosphereModel,
    substepper,
    Δt,
    β_stage,
    Uᴸ,
    advection
)

Execute one Wicker–Skamarock RK3 stage of the linearized acoustic substep loop. Number and size of substeps in this stage depend on substepper.substep_distribution.

freeze_linearization_state!(
    substepper::AcousticSubstepper,
    model
)

Compute the background quantities used by the substepper as the first linearization point of an outer step. Subsequent RK stages call prepare_acoustic_cache!, which refreshes the same cached quantities to the stage-entry state.

After this call:

• linearization_exner = Πᴸ = (pᴸ/pˢᵗ)^κ derived from model.dynamics.pressure
• linearization_potential_temperature = θᴸ = ρθᴸ/ρᴸ derived from model.dynamics.dry_density + ρθ

prepare_acoustic_cache!(
    substepper::AcousticSubstepper,
    model
)

Stage-start cache preparation. Refreshes the cached linearization quantities (Πᴸ, θᴸ, γᵐRᵐᴸ) to the stage-entry state $Uᴸ_\mathrm{stage}$ (per Skamarock & Klemp 2008 above eq. 16), recomputing them from the live model.dynamics.*. The rewind-perturbation initialization (initialize_stage_perturbations!, called next) handles the WS-RK3 invariant by setting $(ρ)′_\mathrm{init} = Uᴸ_\mathrm{outer} − Uᴸ_\mathrm{stage}$ (zero for stage 1; nonzero for stages 2 and 3).

stage_fractions(
    _::WickerSkamarock3
) -> Tuple{Rational{Int64}, Rational{Int64}, Rational{Int64}}

Return the stage-fraction tuple $(β_1, β_2, β_3)$ for the outer Runge–Kutta scheme. For WickerSkamarock3 this is the canonical $(1/3, 1/2, 1)$ of Wicker and Skamarock (2002).

SpecificForcing(
    forcing
) -> SpecificForcing{_A, Nothing, Nothing} where _A

Wrap a user-supplied forcing that produces a specific (per-unit-mass) tendency so that Breeze applies the density multiply $ρ$ at kernel time. After materialization, the kernel callable returns

\[ρ(i, j, k) \, F_ϕ(i, j, k, t)\]

interpolating $ρ$ to the appropriate cell face for fields whose target prognostic lives at Face in any direction (e.g. $ρ$ is interpolated to x-Face for u-forcings via $ℑxᶠᵃᵃ$, to z-Face for w via $ℑzᵃᵃᶠ$). Under AnelasticDynamics, ρ is the reference density ρᵣ(z). Under CompressibleDynamics, the carrier depends on the target conservation law: momentum and thermodynamic tendencies use the dry-air coupling density ρᵈ, while moisture, microphysical moments, and user tracers use total density ρ. The same wrapper handles all carriers.

Users typically supply specific forcings directly through specific-named keys (u, v, w, θ, s, qᵉ, qᵛ, …) in the forcing NamedTuple passed to AtmosphereModel, and the dispatch wraps each entry in SpecificForcing automatically. The wrapper can also be constructed directly when finer control is needed.

The inner forcing can be anything accepted by Breeze's materialize_atmosphere_model_forcing: a function (x, y, z, t), a Returns callable, a Field, an Oceananigans.Forcing, a Breeze forcing such as SubsidenceForcing or one produced by geostrophic_forcings, or a tuple of these.

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

Forcing that represents large-scale subsidence advecting horizontally-averaged fields downward. The kernel returns the specific tendency

\[F_ϕ = - w^s \, ∂_z \overline{ϕ}\]

where $w^s$ is the subsidence_vertical_velocity and $\overline{ϕ}$ is the horizontal average of the field being forced. Supply SubsidenceForcing under the specific prognostic name (e.g. θ, qᵉ, u); the AtmosphereModel dispatch wraps it in SpecificForcing so the density factor $ρ$ is applied automatically at kernel time.

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 Breezegrid = 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.ρθ.forcing# outputSubsidenceForcing 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, keyed under specific names u and v. Each GeostrophicForcing returns a specific tendency; the model's density factor ρ is applied automatically via SpecificForcing when the forcing is dispatched under a specific key, with the correct horizontal interpolation of ρ to the appropriate cell face.

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 Breezeuᵍ(z) = -10 + 0.001zvᵍ(z) = 0.0coriolis = FPlane(f=1e-4)forcing = geostrophic_forcings(uᵍ, vᵍ)# outputNamedTuple with 2 GeostrophicForcings:├── u: GeostrophicForcing{XDirection}│   └── geostrophic_velocity: vᵍ (generic function with 1 method)└── v: GeostrophicForcing{YDirection}    └── geostrophic_velocity: 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, SP, 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,
    base_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. base_pressure is the pressure datum at z = 0. On a raised domain, the default bottom-face pressure extends the lowest prescribed density down to the datum; pass surface_pressure to prescribe a different hydrostatic anchor explicitly.

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

The grid's bottom face is the ground, so give it as z = (0, Lz); on a domain that starts at $z = 0$ the bottom-face surface_pressure and the $z = 0$ base_pressure datum coincide.

using Oceananigansusing Breezegrid = RectilinearGrid(size=(4, 4, 8), x=(0, 1000), y=(0, 1000), z=(0, 2000))reference_state = ReferenceState(grid, ThermodynamicConstants())dynamics = PrescribedDynamics(reference_state)(dynamics.surface_pressure[1, 1, 1], dynamics.base_pressure)# output(101325.0, 101325.0)

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

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
• negative_moisture_correction: Correction scheme for negative moisture produced by advection.
  • nothing (default): No correction
    • VerticalBorrowing(): Vertical redistribution of the moisture prognostic only
  • SpeciesBorrowing(): Same-level species borrowing only
  • SpeciesBorrowing(vertical_borrowing=VerticalBorrowing()): Species borrowing with vertical redistribution

struct BulkMicrophysics{N, C, B, NMC}

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)
• negative_moisture_correction: Correction scheme for negative moisture produced by advection.
  • nothing (default): No correction
    • VerticalBorrowing(): Vertical redistribution of the moisture prognostic only
  • SpeciesBorrowing(): Same-level species borrowing only
  • SpeciesBorrowing(vertical_borrowing=VerticalBorrowing()): Species borrowing with vertical redistribution

struct ConstantRateCondensateFormation{FT} <: Breeze.Microphysics.AbstractCondensateFormation

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.

InstantaneousPrecipitation(
;
    ...
) -> InstantaneousPrecipitation{S} where S<:(SaturationAdjustment{WarmPhaseEquilibrium, S} where S<:SecantSolver)
InstantaneousPrecipitation(
    FT::DataType;
    equilibrium,
    solver
) -> InstantaneousPrecipitation{S} where S<:(SaturationAdjustment{WarmPhaseEquilibrium, S} where S<:SecantSolver)

Construct an InstantaneousPrecipitation scheme. equilibrium selects the phase equilibrium used by the underlying saturation solve (default warm-phase), and solver controls its iteration (see SaturationAdjustment).

struct InstantaneousPrecipitation{S}

Instantaneous-precipitation microphysics: an instantaneous, irreversible condensation with immediate rain-out and no re-evaporation (no cloud or rain stage). Excess water vapor above saturation condenses, releases its latent heat to the air, and is removed as precipitation in the same step. This is the "large-scale condensation" of the DCMIP2016 Reed–Jablonowski simple-physics suite.

The saturation/latent-heat solve is delegated to a SaturationAdjustment instance (saturation_adjustment), so the equilibrium thermodynamics are shared and validated. The distinction is irreversibility: the condensate is purged from the prognostic vapor every step and cannot re-evaporate.

The prognostic moisture is the vapor density ρqᵛ (no condensate is retained).

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, S} where {E<:MixedPhaseEquilibrium, S<:SecantSolver}
SaturationAdjustment(
    FT::DataType;
    solver,
    equilibrium,
    tolerance,
    maxiter
) -> SaturationAdjustment{E, S} where {E<:MixedPhaseEquilibrium, S<:SecantSolver}

Return SaturationAdjustment microphysics representing an instantaneous adjustment to equilibrium between condensates and water vapor, computed by a secant iteration on the temperature residual controlled by solver.

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.

The options for solver are SecantSolver (default: SecantSolver(abstol=1e-4, maxiter=20), an absolute tolerance on the temperature-like residual in Kelvin) and FixedIterations, which performs a fixed number of secant steps with no convergence test (the form required for Reactant tracing and cheap reverse-mode differentiation).

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

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 Breezegrid = RectilinearGrid(size=(1, 1, 128), extent=(1e3, 1e3, 1e3))microphysics = SaturationAdjustment()model = AtmosphereModel(grid; microphysics)set!(model, θ=300, qᵗ=0.005)  # subsaturatedℋ = RelativeHumidity(model)# outputKernelFunctionOperation 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# output1×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.214949, min=0.137169, mean=0.172626

We also provide a convenience constructor for the Field:

ℋ_field = RelativeHumidityField(model)# output1×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.214949, min=0.137169, mean=0.172626

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

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

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

Saturation adjustment for the LiquidIceDensityState: a secant on the constant-density θˡⁱ-conservation residual, so qsat and the θˡⁱ inversion are evaluated at the state's actual density ρ (with true pressure p = ρRᵐT) rather than a fixed reference pressure. This is the density-consistent analogue of the generic (reference-pressure) adjust_state secant; like that one it holds θˡⁱ fixed (conserves it). See NumericalEarth/Breeze.jl#765.

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 the terminal_velocity_coefficient, $Cᵨ$ is the density_scale, and $β^𝕎$ is the terminal_velocity_exponent.

number_concentration(model, species::Symbol) -> Any

Lazy diagnostic returning the total number concentration $ρnˣ$ (m⁻³) for the requested species.

For OneMomentCloudMicrophysics, species ∈ (:rain, :snow) returns a KernelFunctionOperation that computes $n_0 \, λ^{-1}$ from the prognostic $ρqˣ$ and the scheme's size distribution. Snow's intercept $n_0$ depends on $(q, ρ)$ per Kaul et al. (2015) — so this diagnostic stays consistent with the scheme's actual DSD without re-encoding species-specific physics at every call site.

For TwoMomentCloudMicrophysics, returns the prognostic $ρnˣ$ field directly (e.g., :rain → ρnʳ, :cloud_liquid → ρnᶜˡ).

For PredictedParticlePropertiesMicrophysics, returns the prognostic $ρnˣ$ field when present. In the default prescribed-cloud-number configuration, :cloud_liquid returns a lazy constant operation based on the configured droplet number.

Returns nothing if the species is not carried by the model (e.g., :hail for a 1-mom scheme without hail). Errors for microphysics schemes that do not define a DSD-based number concentration (e.g., SaturationAdjustment).

The return shape is therefore polymorphic: a lazy KernelFunctionOperation for diagnosed or prescribed concentrations and a stored Field for prognostic concentrations. Use number_concentration_field when you want a uniformly Field-typed handle.

number_concentration_field(model, species::Symbol) -> Any

Field-typed handle for the number_concentration diagnostic. Allocates a Field shell around lazy KernelFunctionOperations (use compute! to populate it) and returns prognostic $ρnˣ$ fields directly. Returns nothing when the requested species is not carried by the model.

PredictedParticleProperties

Predicted Particle Properties (P3) microphysics scheme implementation.

P3 is a bulk microphysics scheme that uses a single ice category with continuously predicted properties (rime fraction, rime density, liquid fraction) rather than multiple discrete ice categories.

Key Features

• Single ice category with predicted properties
• Two-moment ice (mass and number)
• Predicted liquid fraction on ice particles
• Rime fraction and rime density evolution
• Ice-side integrals are read from the ASCII lookup tables; rain 1D integrals are tabulated at startup using Chebyshev–Gauss quadrature

Complete Reference List

This implementation is based on the following P3 papers:

1. Morrison & Milbrandt (2015a) - Original P3: $m(D)$, $A(D)$, $\mathbb{W}(D)$, process rates Morrison and Milbrandt (2015a)

2. Morrison et al. (2015b) - Part II: Case study validation Morrison et al. (2015b)

3. Milbrandt & Morrison (2016) - Part III: Multiple ice categories (NOT implemented) Milbrandt and Morrison (2016)

4. Milbrandt et al. (2025) - Predicted liquid fraction: shedding, refreezing Milbrandt et al. (2025)

Source Code

Based on P3-microphysics v5.5.0

Not Implemented

• Three-moment ice (prognostic reflectivity) from Milbrandt et al. (2021)
• Multiple free ice categories from Milbrandt & Morrison (2016)
• Lookup table I/O for all table types

AerosolActivation(
    mode1::Breeze.Microphysics.PredictedParticleProperties.AerosolMode{FT},
    rest::Breeze.Microphysics.PredictedParticleProperties.AerosolMode{FT}...;
    prognostic,
    thermodynamic_constants,
    molecular_weight_water,
    universal_gas_constant,
    activation_timescale,
    surface_tension_reference,
    surface_tension_temperature_derivative,
    surface_tension_reference_temperature,
    lognormal_activation_factor,
    activated_droplet_radius,
    activation_supersaturation_threshold,
    minimum_supersaturation,
    minimum_saturation_mass_fraction
) -> Breeze.Microphysics.PredictedParticleProperties.AerosolActivation{_A, _B, <:Tuple{Breeze.Microphysics.PredictedParticleProperties.AerosolMode, Vararg{Breeze.Microphysics.PredictedParticleProperties.AerosolMode}}} where {_A, _B}

Construct an AerosolActivation from one or more AerosolModes.

The activation timescale $τ_{act}$ controls how quickly the cloud droplet number relaxes toward the activated equilibrium. Default 1.0 s.

Passing an AerosolActivation to P3 makes cloud droplet number prognostic. The prognostic keyword is separate: true carries the unactivated aerosol reservoir $ρn^a$, which activation depletes, while the default false holds the population fixed at the distribution total.

Everything else the activation physics needs is a keyword here rather than a literal in activated_number: the condensate density and molecular weight, the linear surface-tension fit $σ(T) = σ_0 + (dσ/dT)(T - T_{ref})$ whose defaults are water's, the $3\sqrt{2}$ factor of the lognormal activation integral (lognormal_activation_factor, held at the conventional rounded 4.242), the radius of a newly activated droplet, the supersaturation above which activation proceeds, and two floors that keep the supersaturation finite. A condensable species other than water is configured by overriding the first group here and in AerosolMode.

By default, the water molecular weight, liquid-water density, universal gas constant, and surface-tension reference temperature come from thermodynamic_constants.

Examples

using Breeze.Microphysics.PredictedParticleProperties: AerosolActivation, AerosolModeaerosol = AerosolActivation(AerosolMode())length(aerosol.modes)# output1

using Breeze.Microphysics.PredictedParticleProperties: AerosolActivation, AerosolModeaerosol = AerosolActivation(    AerosolMode(number_mixing_ratio=100e6, mean_radius=0.08e-6),    AerosolMode(number_mixing_ratio=50e6,  mean_radius=1.0e-6, geometric_std=2.5);    activation_timescale = 2.0)length(aerosol.modes)# output2

The aerosol population is fixed by default:

using Breeze.Microphysics.PredictedParticleProperties: AerosolActivation, AerosolModesummary(AerosolActivation(AerosolMode()))# output"AerosolActivation(1 mode, fixed reservoir)"

AerosolMode(
;
    ...
) -> Breeze.Microphysics.PredictedParticleProperties.AerosolMode
AerosolMode(
    FT::DataType;
    number_mixing_ratio,
    mean_radius,
    geometric_std,
    vant_hoff_factor,
    osmotic_potential,
    mass_fraction_soluble,
    aerosol_density,
    molecular_weight_aerosol,
    thermodynamic_constants,
    molecular_weight_water
) -> Breeze.Microphysics.PredictedParticleProperties.AerosolMode

Construct an AerosolMode representing one component of a multimodal aerosol size distribution. Particles in a mode share one chemical composition and their radii follow a lognormal distribution described by mean_radius and geometric_std. Multiple modes can therefore represent distinct aerosol populations, such as Aitken and accumulation particles.

The solute activity parameter $β_{act} = ν_i ϕ_s ε_m M_w ρ_a / (M_a ρ_w)$ is precomputed at construction time from the chemistry parameters.

Default chemistry is ammonium sulfate (NH₄)₂SO₄.

Keyword Arguments

• number_mixing_ratio: Aerosol number per unit mass of air [kg⁻¹], default 300×10⁶. With prognostic=true in AerosolActivation, the reservoir ρnᵃ is initialized to air density times the total over all modes once density is available. set! repeats this initialization unless nᵃ [kg⁻¹] or ρnᵃ [m⁻³] is supplied.
• mean_radius: Geometric mean radius [m], default 0.05 μm
• geometric_std: Geometric standard deviation [-], default 2
• vant_hoff_factor: van't Hoff factor [-], default 3
• osmotic_potential: Osmotic potential [-], default 1
• mass_fraction_soluble: Mass fraction soluble [-], default 0.9
• aerosol_density: Aerosol density [kg/m³], default 1777
• molecular_weight_aerosol: Molecular weight of aerosol [kg/mol], default 0.132
• thermodynamic_constants: Constants supplying the water molecular weight and liquid-water density
• molecular_weight_water: Molecular weight of the condensate [kg/mol], default thermodynamic_constants.vapor.molar_mass

The molecular weight and thermodynamic_constants.liquid.density enter only through $β_{act}$. A condensable species other than water is configured through thermodynamic_constants and the surface-tension fit in AerosolActivation.

References

Morrison and Grabowski (2007)

Examples

using Breeze.Microphysics.PredictedParticleProperties: AerosolModemode = AerosolMode()mode.mean_radius# output5.0e-8

CloudDroplets(
;
    ...
) -> Breeze.Microphysics.PredictedParticleProperties.CloudDroplets
CloudDroplets(
    FT;
    number_concentration,
    condensation_timescale,
    shape,
    shape_parameter
) -> Breeze.Microphysics.PredictedParticleProperties.CloudDroplets

Construct CloudDroplets with prescribed parameters.

Cloud droplets in P3 are treated simply: their number concentration is prescribed rather than predicted. This is a common simplification appropriate for many applications where aerosol-cloud interactions are not the focus.

Why prescribe Nᶜˡ?

Predicting cloud droplet number Nᶜˡ requires treating aerosol activation physics, which adds substantial complexity. For simulations focused on ice processes or bulk precipitation, prescribed Nᶜˡ is sufficient.

The prescribed-Nᶜˡ simplification means: (1) homogeneous freezing below −40°C transfers the prescribed Nᶜˡ rather than a locally depleted droplet count, and (2) autoconversion sensitivity to Nᶜˡ is controlled by the prescribed value rather than dynamically. Pass aerosol = AerosolActivation(AerosolMode()) to predict Nᶜˡ instead.

There is no separate mass-number consistency cap on homogeneous freezing: homogeneous_freezing_cloud_rate transfers all of Nᶜˡ at T < homogeneous_freezing_temperature. compute_p3_process_rates diagnoses the rate from the post-process residual cloud and rescales mass and number together by a single sink_limiting_factor. This limits the frozen mass to the residual cloud while preserving its diagnosed mass-number ratio; it does not impose a minimum frozen-particle mass or independently limit the transferred number.

Cloud DSD shape parameter: Process rates diagnose μᶜˡ from the local droplet concentration via liu_daum_shape_parameter, using the relation and bounds in shape. This applies to both prescribed and prognostic droplet number. Configure shape for sensitivity studies that change the simulated PSD.

The shape_parameter keyword replaces the value stored in the shape_parameter field. It does not change the PSD or freezing rates used in a simulation.

The stored freezing_psd_correction = Γ(μᶜˡ+7)Γ(μᶜˡ+1)/Γ(μᶜˡ+4)² is evaluated at shape_parameter. immersion_freezing_cloud_rate recomputes the correction from the local μᶜˡ with psd_correction_spherical_volume.

Typical values:

• Continental: Nᶜˡ ~ 100-300 × 10⁶ m⁻³ → μᶜˡ ~ 4–8
• Marine: Nᶜˡ ~ 50-100 × 10⁶ m⁻³ → μᶜˡ ~ 8–10

Autoconversion: Cloud droplets are converted to rain via collision-coalescence following Khairoutdinov and Kogan (2000).

Keyword Arguments

• number_concentration: Nᶜˡ [1/m³], default 200×10⁶
• condensation_timescale: Saturation relaxation [s], default 1.0. The coupled adjustment uses ProcessRate.sink_limiting_timescale instead.
• shape: CloudShape holding the coefficients and bounds of the Liu-Daum relation, default CloudShape(FT). Read by every path that diagnoses μᶜˡ from a local droplet number.
• shape_parameter: μᶜˡ for cloud gamma PSD [-], default nothing (diagnosed from Nᶜˡ via Liu-Daum relation). Pass an explicit value to override the construction-time diagnosis only.

References

Morrison and Milbrandt (2015a), Khairoutdinov and Kogan (2000).

Examples

using Oceananigans, Breezeusing Breeze.Microphysics.PredictedParticleProperties: CloudDropletscloud = CloudDroplets()round(cloud.shape_parameter, digits=1)  # μᶜˡ diagnosed from Nᶜˡ = 200×10⁶ m⁻³# output5.7

CloudShape(
;
    ...
) -> Breeze.Microphysics.PredictedParticleProperties.CloudShape
CloudShape(
    FT::DataType;
    relative_dispersion_number_coefficient,
    relative_dispersion_intercept,
    minimum_shape_parameter,
    maximum_shape_parameter
) -> Breeze.Microphysics.PredictedParticleProperties.CloudShape

Construct CloudShape.

Keyword Arguments

• relative_dispersion_number_coefficient: $\mathbb{C}_{cl,1}$ [m³], default 5.714e-10
• relative_dispersion_intercept: $\mathbb{C}_{cl,2}$ [-], default 0.2714
• minimum_shape_parameter: $\mathbb{C}_{cl,3}$ [-], default 2
• maximum_shape_parameter: $\mathbb{C}_{cl,4}$ [-], default 15

Examples

using Breeze.Microphysics.PredictedParticleProperties: CloudShapeCloudShape(Float64)# outputCloudShape(ℂᶜˡ₁=5.714e-10 m³, ℂᶜˡ₂=0.2714, ℂᶜˡ₃=2.0, ℂᶜˡ₄=15.0)

CloudShape{FT}

Coefficients of the Liu-Daum (2000)-type relation that diagnoses the cloud gamma PSD shape parameter $μ^{cl}$ from the droplet number density, together with the bounds the diagnosis is clamped to. Evaluated by liu_daum_shape_parameter:

\[\chi = \mathbb{C}_{cl,1} \, N^{cl} + \mathbb{C}_{cl,2}, \qquad \mu^{cl} = \mathrm{clamp}\!\left(\frac{1}{\chi^2} - 1,\; \mathbb{C}_{cl,3},\; \mathbb{C}_{cl,4}\right)\]

$\chi$ is the relative dispersion of the droplet spectrum, and $(\mathbb{C}_{cl,1}, \mathbb{C}_{cl,2})$ are the Liu-Daum regression of $\chi$ on droplet concentration, fit to aircraft measurements of warm cloud droplet spectra: at fixed water content, more droplets means a narrower spectrum, hence a larger $μ^{cl}$. The free parameters $(\mathbb{C}_{cl,3}, \mathbb{C}_{cl,4})$ bound $μ^{cl}$ to the range over which the fit was measured.

The coefficient is stated here for the absolute number density in SI units [m⁻³], so $\mathbb{C}_{cl,1}$ carries units of m³. The published form uses cm⁻³, hence the 10⁻⁶ difference from the printed 5.714 × 10⁻⁴.

See the constructor for the meaning, units and defaults of each coefficient.

IceBulk(
;
    ...
) -> Breeze.Microphysics.PredictedParticleProperties.IceBulk{_A, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing} where _A
IceBulk(
    FT::DataType;
    maximum_mean_diameter,
    minimum_mean_diameter
) -> Breeze.Microphysics.PredictedParticleProperties.IceBulk{_A, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing} where _A

Construct IceBulk with parameters and quadrature-based integrals.

These integrals compute bulk properties by averaging over the particle size distribution. They are used for radiation, radar, and diagnostics.

Diagnostic integrals:

• effective_radius: Radiation-weighted radius $r_e = (3/(4ρ_i^*)) ∫m·N'dD / ∫A·N'dD$, with the table generator's reference ice density $ρ_i^* = 916.7$ kg/m³
• mean_diameter: Mass-weighted diameter $D_m = ∫D·m·N'dD / ∫m·N'dD$
• mean_density: Mass-weighted density $ρ̄ = ∫ρ·m·N'dD / ∫m·N'dD$
• reflectivity: Number-normalized equivalent radar reflectivity. Dry ice uses $0.1892 ∫D_{eq}^6 N'dD / ∫N'dD$, with equivalent diameter computed from particle mass at 917 kg/m³. Partially melted ice uses the generator's wet-ice scattering calculation; the fully liquid limit uses $D^6$. Multiply by ice number density for the volume integral.

Tabulated distribution parameters:

• slope: Slope parameter λ recorded when the table was generated
• shape: Shape parameter μⁱ recorded when the table was generated, read back by compute_ice_shape_parameter

Process integrals:

• shedding: Rate at which meltwater sheds from large particles

Keyword Arguments

• maximum_mean_diameter: Upper Dm limit [m], default 0.02 (2 cm)
• minimum_mean_diameter: Lower Dm limit [m], default 2×10⁻⁶ (2 μm)

References

Morrison and Milbrandt (2015a), Field et al. (2007) for μⁱ-λ relationship.

IceCollection(

) -> Breeze.Microphysics.PredictedParticleProperties.IceCollection{Nothing, Nothing, Nothing, Nothing}

Construct IceCollection with placeholder (nothing) integrals, following the materialization pattern: read_lookup_tables replaces each field with the corresponding tabulated integral.

Collection processes describe ice particles sweeping up other hydrometeors through gravitational settling. The integrals held here are the ones that depend on the ice size distribution alone:

Aggregation (ice + ice → larger ice): Ice particles collide and stick together to form larger aggregates. This is the dominant growth mechanism for snow, and depends on the differential fall speeds of particles of different sizes. Consumed by ice_aggregation_rate.

Cloud collection (ice + cloud droplets → rime on ice): The PSD-integrated sweep-out kernel $\int \mathbb{W}(D) A(D) N'(D) \, dD$ [m³/s] per particle, with the collision kernel set to zero for ice diameters below 100 μm. Cloud droplets are small enough relative to ice that their own size distribution does not enter the collision geometry, so a single ice-PSD integral suffices. Consumed by cloud_riming_rate below freezing and by cloud_warm_collection_rate above it.

Aerosol scavenging (cloud_aerosol_collection, ice_aerosol_collection): Collection by ice particles of water-friendly and ice-friendly interstitial aerosol, respectively.

Ice-rain collection is handled separately, by IceRainCollection and the 5D rain-ice block embedded in Lookup Table 1, because its kernel needs the rain slope parameter $λ_r$ in addition to the ice PSD.

Collection efficiencies are not stored here. They live in ProcessRate alongside the other rate parameters, as cloud_ice_collection_efficiency ($E^{ci}$) and rain_ice_collection_efficiency ($E^{ri}$).

References

Morrison and Milbrandt (2015a) Sections 2d-e, Milbrandt and Yau (2005).

IceDeposition(

) -> Breeze.Microphysics.PredictedParticleProperties.IceDeposition{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}
IceDeposition(
    ::DataType
) -> Breeze.Microphysics.PredictedParticleProperties.IceDeposition{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}

Construct IceDeposition with quadrature-based ventilation integrals.

Ice growth/decay by vapor deposition/sublimation follows the diffusion equation with ventilation enhancement. The ventilation factor $fᵛᵉ$ accounts for enhanced vapor transport due to particle motion through air:

\[fᵛᵉ = a + b \cdot Sc^{1/3} Re^{1/2}\]

where $Sc$ is the Schmidt number and $Re$ is the Reynolds number. Hall and Pruppacher (1976) showed that falling particles have significantly enhanced vapor exchange compared to stationary particles.

Thermal conductivity $κ$ and vapor diffusivity $Dᵥ$ are computed at runtime from temperature, pressure, and the model thermodynamic constants via air_transport_properties(T, P, constants). They are not stored on IceDeposition.

Basic ventilation integrals:

• ventilation: Integrated over full size spectrum
• enhanced_ventilation: For larger particles (D > 100 μm)

Size-regime ventilation (for melting with liquid fraction):

• small_ice_ventilation_*: D ≤ Dcrit, meltwater → rain
• large_ice_ventilation_*: D > Dcrit, meltwater → liquid on ice

References

Hall and Pruppacher (1976), Morrison and Milbrandt (2015a) Eq. 34.

IceFallSpeed(
;
    ...
) -> Breeze.Microphysics.PredictedParticleProperties.IceFallSpeed{_A, Nothing, Nothing} where _A
IceFallSpeed(
    FT::DataType;
    thermodynamic_constants,
    reference_pressure,
    reference_temperature,
    reference_air_density
) -> Breeze.Microphysics.PredictedParticleProperties.IceFallSpeed{_A, Nothing, Nothing} where _A