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.

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

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 surface_pressure (a scalar; defaults to the dynamics' mean surface pressure). For CompressibleDynamics.

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.

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

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! invokes microphysics. 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)$. For compressible dynamics, this is the prognostic pressure field. The anomaly and total-pressure counterparts are pressure_anomaly and total_pressure.

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 total initial aerosol number concentration [m⁻³] for a microphysics scheme.

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

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

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, surface_pressure, thermodynamic_constants,
                                                microphysical_fields, specific_prognostic_moisture, temperature)

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

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

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

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

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

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

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

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

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

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. The recomputed reference is horizontally uniform (a single column). (Assumes a 1-D column reference, the form built from a constant or z-dependent θ₀.)

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 (ρe, ρqᵗ, ρu, etc.) are rescaled by ρᵣ_new / ρᵣ_old so that the specific quantities (e, 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.

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 (ρθ, ρe, ρ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., :ρθ, :ρe, :ρE). Accepts a Symbol, Val(Symbol), or formulation struct.

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

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

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

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

This is the prognostic thermodynamic variable used in StaticEnergyThermodynamics.

Arguments

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

Examples

using Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300)e = StaticEnergy(model)Field(e)# 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 surface 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 (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ᴰ.

When a FilteredSurfaceVelocities is supplied via filtered_velocities, every field entering the formula — the wind speed |U|, the velocity u, and the virtual potential temperature θᵥ 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)).

Monin–Obukhov consistency

ρ₀ is computed from surface quantities (surface_pressure and surface_temperature) via the ideal gas law, so it is a true surface density — independent of the vertical grid resolution. Using the prognostic density at the first cell would introduce a grid-dependent ρ₀ (the first-cell height ½Δz shifts the value as the grid is refined), which is inconsistent with the bulk-flux closure derived from Monin–Obukhov similarity.

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() or YDirection()). If nothing, the direction is inferred from the field location during boundary condition regularization.
• coefficient: The drag coefficient (default: 1e-3). Can be a constant or a PolynomialCoefficient for wind and stability-dependent transfer coefficients.
• gustiness: Minimum wind speed to prevent singularities when winds are calm (default: 0)
• surface_temperature: Surface temperature, used to compute ρ₀ and required when using PolynomialCoefficient with stability correction. Can be a Field, Function, or Number. (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, and $Δϕ$ is the difference between the near-surface atmospheric value and the surface value of the thermodynamic variable appropriate to the formulation:

• For LiquidIcePotentialTemperatureFormulation: $Δϕ = θ - θ₀$, where $θ₀ = T₀ / Π₀$ and $Π₀ = (p₀ / pˢᵗ)^{Rᵈ / cᵖᵈ}$ (potential temperature flux)
• For StaticEnergyFormulation: $Δϕ = e - cᵖᵈ T₀$ (static energy flux)

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

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

Keyword Arguments

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

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

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

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

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

Keyword Arguments

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

EnergyFluxBoundaryConditionFunction

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

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

The relationship is:

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

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

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

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 virtual potential temperature 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 θ̄ᵥ 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.

θ̄ᵥ is updated from the virtual-potential-temperature diagnostic owned by any attached PolynomialCoefficient; when the bulk coefficient is a plain Number (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 read 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_parameters = StabilityFunctionParameters(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_parameters: StabilityFunctionParameters (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.StabilityFunctionParameters})}
PolynomialCoefficient(
    FT;
    polynomial,
    roughness_length,
    minimum_wind_speed,
    stability_function,
    surface,
    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.StabilityFunctionParameters})}

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) × 10^{-3}\]

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

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

When polynomial is nothing, the appropriate Large 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) for momentum
• BulkSensibleHeatFlux: default_neutral_sensible_heat_polynomial = (0.128, 0.068, 2.43) for sensible heat
• BulkVaporFlux: default_neutral_latent_heat_polynomial = (0.120, 0.070, 2.55) for latent heat

Keyword Arguments

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

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

Examples

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

using Breeze.BoundaryConditions: PolynomialCoefficient# With explicit polynomialcoef = PolynomialCoefficient(polynomial = (0.142, 0.076, 2.7))# outputPolynomialCoefficient{Float64}├── polynomial: (0.142, 0.076, 2.7)├── roughness_length: 0.00015 m├── minimum_wind_speed: 0.1 m/s├── surface: PlanarLiquidSurface└── stability_function: FittedStabilityFunction (Li et al. 2010)

using Breeze.BoundaryConditions: PolynomialCoefficient# No stability correctioncoef = PolynomialCoefficient(stability_function = nothing)# outputPolynomialCoefficient{Float64}├── polynomial: nothing├── roughness_length: 0.00015 m├── minimum_wind_speed: 0.1 m/s├── surface: PlanarLiquidSurface└── 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.

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

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

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

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

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

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

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

References

• 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ᵖᵐ$.

The relationship is:

\[𝒬 = Jᶿ cᵖᵐ\]

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

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

Create a FluxBoundaryCondition for surface momentum drag.

See BulkDragFunction for details.

Examples

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

BulkSensibleHeatFlux(; coefficient, gustiness=0, surface_temperature)

Create a FluxBoundaryCondition for surface sensible heat flux.

The bulk formula computes

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

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

See BulkSensibleHeatFluxFunction for details.

Example

using BreezeT₀(x, y) = 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, gustiness=0)

Create a FluxBoundaryCondition for surface moisture flux.

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

See BulkVaporFluxFunction for details.

Example

using BreezeT₀(x, y) = 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ᵖᵐ$ to obtain the potential temperature flux: $Jᶿ = 𝒬 / cᵖᵐ$.

ThetaFluxBoundaryCondition(flux)

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

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

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

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

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

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

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

The solar zenith angle $θ_z$ satisfies:

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

where:

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

Arguments

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

Returns

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

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

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

equation_of_time(day_of_year) -> Any

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

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

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

References

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

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

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

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

Arguments

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

solar_declination(day_of_year) -> Any

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

Uses the approximation by Spencer (1971):

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

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

References

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

CompressibleEquations

Module implementing fully compressible dynamics for atmosphere models.

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

The fully compressible Euler equations in conservation form are:

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

Pressure is computed from the ideal gas law:

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

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

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}

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: BatchedTridiagonalSolver for the implicit (ρw)′ update.

AcousticSubstepper(
    grid,
    split_explicit::SplitExplicitTimeDiscretization;
    prognostic_momentum,
    substep_floattype
) -> AcousticSubstepper{_A, _B, _C, AD, _D, CF, MP, TAV, GT, TS} 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})}

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)
• surface_pressure: Mean pressure at the bottom of the atmosphere p₀
• 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,
    surface_pressure,
    reference_potential_temperature,
    reference_temperature,
    reference_vapor_mass_fraction,
    slope_stencil,
    terrain_metrics,
    reference_state,
    temperature_tolerance,
    temperature_maxiter
) -> 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)

• surface_pressure: Mean surface pressure (default: 101325.0 Pa)

• reference_potential_temperature: Potential temperature for building a fixed hydrostatically-balanced reference state used in base-state subtraction. Can be a constant θ₀ or a function θ(z). Default: nothing, 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ᴸ
)

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ᵃᵃᶠ$). ρ is ρᵣ(z) under AnelasticDynamics and the prognostic ρ(x, y, z, t) under CompressibleDynamics; the same wrapper handles both.

Users typically supply specific forcings directly through specific-named keys (u, v, w, θ, e, 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, FT}

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

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

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

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

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

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

Example

using Oceananigansusing Breezegrid = RectilinearGrid(size=(4, 4, 8), extent=(1000, 1000, 2000))reference_state = ReferenceState(grid, ThermodynamicConstants())dynamics = PrescribedDynamics(reference_state)# outputPrescribedDynamics├── density: PrescribedDensity├── pressure: 1×1×8 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on RectilinearGrid on CPU├── surface_pressure: 101325.0└── standard_pressure: 100000.0

BulkMicrophysics(
;
    ...
) -> BulkMicrophysics{N, Nothing, Nothing, 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

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ᶜˡ).

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 1-mom and a stored Field for 2-mom — so the function is snake-cased rather than PascalCased. 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. For 1-mom, allocates a Field shell around the lazy KernelFunctionOperation (use compute! to populate it). For 2-mom, returns the prognostic $ρnˣ$ field directly. Returns nothing when the requested species is not carried by the model.

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

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

Example

using Breeze, Oceananigansgrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 3e3))buoyancy = MoistAirBuoyancy(grid)# outputMoistAirBuoyancy:├── reference_state: ReferenceState{Float64}(p₀=101325.0, θ₀=288.0, pˢᵗ=100000.0)└── thermodynamic_constants: ThermodynamicConstants{Float64}

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

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

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

Construct ParcelDynamics with default (uninitialized) state.

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

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

Lagrangian parcel dynamics for AtmosphereModel.

Fields

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

mutable struct ParcelInitialState{FT, MP}

• x::Any

• y::Any

• z::Any

• w::Any

• qᵗ::Any

• ℰ::Any

• μ::Any

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

ParcelModel

Type alias for AtmosphereModel{<:ParcelDynamics}.

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

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

See also ParcelDynamics, AtmosphereModel.

mutable struct ParcelState{FT, TH, MP}

• x::Any

• y::Any

• z::Any

• w::Any

• ρ::Any

• qᵗ::Any

• ρqᵗ::Any

• ℰ::Any

• ρℰ::Any

• 𝒰::Any

• μ::Any

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

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

• w: parcel vertical velocity [m/s], prognostic for PrognosticVerticalVelocity, zero for PrescribedVerticalVelocity
• ρ: environmental density at parcel height [kg/m³], interpolated from background profile, not the parcel's own density. The parcel density is computed from density(𝒰, constants) using the ideal gas law applied to the parcel's thermodynamic state.

mutable struct ParcelTendencies{FT, GM}

• Gx::Any

• Gy::Any

• Gz::Any

• Gw::Any

• Ge::Any

• Gqᵗ::Any

• Gμ::Any

Tendencies (time derivatives) for parcel prognostic variables.

struct ParcelTimestepper{GT, U0, FT}

SSP RK3 time-stepper for ParcelModel.

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

Fields

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

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

Construct a ParcelTimestepper for SSP RK3 time-stepping.

PrescribedVerticalVelocity

Singleton type for prescribed vertical velocity dynamics. The parcel moves following the prescribed environmental vertical velocity field w(z).

This is the default vertical velocity formulation: dz/dt = w_env(z).

PrognosticVerticalVelocity

Singleton type for prognostic vertical velocity dynamics. The parcel has a prognostic vertical velocity driven by buoyancy, i.e., dz/dt = w and dw/dt = b, where b = -g (ρᵖ - ρᵉ) / ρᵉ is the net buoyancy from the density difference, including both the virtual temperature effect and condensate loading.

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

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

Compute tendencies for the parcel prognostic variables.

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

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

compute_vertical_velocity_tendencies!(
    tendencies,
    state,
    dynamics,
    model,
    _::PrescribedVerticalVelocity
)

Compute vertical velocity tendencies for PrescribedVerticalVelocity.

The parcel follows the environmental vertical velocity: dz/dt = w_env(z). The parcel velocity tendency Gw is zero (unused prognostic).

compute_vertical_velocity_tendencies!(
    tendencies,
    state,
    dynamics,
    model,
    _::PrognosticVerticalVelocity
)

Compute vertical velocity tendencies for PrognosticVerticalVelocity.

The parcel has a prognostic vertical velocity driven by buoyancy, i.e., dz/dt = w and dw/dt = B.

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

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

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

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

parcel_buoyancy(state, dynamics, constants) -> Any

Compute the net buoyancy acceleration for a parcel.

The buoyancy is computed from the density difference between the parcel and environment:

\[B = -g (ρ_{parcel} - ρ_{env}) / ρ_{env}\]

Here, $ρ_{env}$ is the environmental density interpolated at the parcel height and $ρ_{parcel} = p / (Rᵐ T)$ is the total parcel density from the ideal gas law, where $Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ$ with $qᵈ = 1 - qᵛ - qˡ - qⁱ$. This formulation captures both the virtual temperature effect (from vapor content) and the condensate loading effect (condensate reduces $qᵈ$, reducing $Rᵐ$, increasing $ρ_{parcel}$) in a single term without double-counting.

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

Apply an SSP RK3 substep with coefficient $α$:

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

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

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

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

Step the parcel state forward using Forward Euler:

\[x^{n+1} = x^n + Δt \, G^n\]

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

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

Copy current prognostic state values to the initial state storage.

PotentialTemperatureFormulations

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

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

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

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

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

Recovering temperature from θˡⁱ may require an iterative inversion, depending on the dynamics: with prognostic-density (compressible) dynamics, temperature solves the implicit relation T = (ρRᵐT/pˢᵗ)^κ θ + ΔL/cᵖᵐ. The inversion is controlled by temperature_solver:

• DefaultTemperatureSolver() (default): resolved at materialization to default_temperature_solver(dynamics) — nothing for anelastic dynamics (closed-form inversion) and NewtonSolver() for compressible dynamics.
• NewtonSolver: tolerance-based Newton iteration.
• FixedIterations: a fixed number of Newton steps with no convergence test, which unrolls to straight-line code (required for Reactant tracing and cheap reverse-mode differentiation).
• nothing: the non-iterated closed-form inversion.

LiquidIcePotentialTemperatureFormulation(
;
    temperature_solver
) -> LiquidIcePotentialTemperatureFormulation{Nothing, Nothing, Breeze.AtmosphereModels.DefaultTemperatureSolver}

Return a LiquidIcePotentialTemperatureFormulation with the given temperature_solver. The prognostic and diagnostic fields are materialized later in the model constructor.

using BreezeLiquidIcePotentialTemperatureFormulation(temperature_solver = FixedIterations(2))# outputLiquidIcePotentialTemperatureFormulation└── temperature_solver: FixedIterations(2)

Solvers

Iterative solvers for the small nonlinear scalar problems that arise in Breeze's thermodynamics: equation-of-state temperature inversions, saturation adjustment, and dewpoint computation.

A "solver" is a lightweight, isbits description of an iteration's stopping rule that algorithms dispatch on:

• NewtonSolver and SecantSolver iterate until a tolerance-based convergence criterion is met (or maxiter is reached).
• FixedIterations performs an exact number of iterations with no convergence test at all. Because the trip count is fixed, the loop unrolls to straight-line code, which is required for Reactant tracing and cheap reverse-mode differentiation (a tolerance-based while loop traces to an XLA while op whose adjoint is pathological — see NumericalEarth/Breeze.jl#767).
• nothing means "do not iterate": the algorithm returns its initial guess (typically a closed-form approximation).

The drivers newton_solve and secant_solve implement the iterations once, so every algorithm shares the same loop logic and the same solver vocabulary.

Tolerance conventions

The choice between reltol and abstol follows the natural scale of the residual:

• Quantities bounded away from zero with a fixed precision target use an absolute tolerance. Every iteration on a temperature (the θˡⁱ→T inversion, saturation adjustment, the Boussinesq adjustment temperature) is solved to abstol = 1e-4 K — far below any physically or numerically meaningful temperature increment, yet reached in a handful of iterations by both Newton (quadratic) and secant (superlinear).
• Quantities that range over orders of magnitude use a relative tolerance against an algorithm-supplied scale. The dewpoint solve iterates on a saturation-vapor-pressure residual (Pa), which spans roughly two decades over the atmospheric temperature range, so it uses reltol = 1e-4 against the vapor pressure; an absolute Pa tolerance would be meaningless across that range.

Iteration caps reflect each method's convergence order and role: the quadratically-convergent Newton inversion caps at maxiter = 8, the superlinear secant temperature solves cap at maxiter = 20, and the dewpoint diagnostic caps at maxiter = 10.

struct FixedIterations

A solver that performs exactly iterations iterations with no convergence test.

The fixed trip count means the iteration unrolls to straight-line, branch-free code, making it the right choice for Reactant tracing and reverse-mode differentiation, where data-dependent while loops are pathological (NumericalEarth/Breeze.jl#767).

using BreezeFixedIterations(2)# outputFixedIterations(2)

NewtonSolver(; ...) -> NewtonSolver
NewtonSolver(
    FT::DataType;
    reltol,
    abstol,
    maxiter
) -> NewtonSolver

Return a NewtonSolver with relative tolerance reltol, absolute tolerance abstol, and iteration cap maxiter.

using BreezeNewtonSolver(maxiter=4)# outputNewtonSolver(reltol=0.0, abstol=0.0001, maxiter=4)

struct NewtonSolver{FT}

A Newton iteration that terminates when the step size Δx satisfies |Δx| ≤ max(abstol, reltol * |x|), or after maxiter iterations.

SecantSolver(; ...) -> SecantSolver
SecantSolver(
    FT::DataType;
    reltol,
    abstol,
    maxiter
) -> SecantSolver

Return a SecantSolver with relative tolerance reltol, absolute tolerance abstol, and iteration cap maxiter.

using BreezeSecantSolver(abstol=1e-4, maxiter=20)# outputSecantSolver(reltol=0.0, abstol=0.0001, maxiter=20)

struct SecantSolver{FT}

A secant iteration that terminates when the residual r satisfies |r| ≤ max(abstol, reltol * |scale|) — where scale is an algorithm-supplied magnitude for the residual — or after maxiter iterations.

StaticEnergyFormulations

Submodule defining the static energy thermodynamic formulation for atmosphere models.

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

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

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

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

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

TerrainFollowingDiscretization

Module implementing terrain-following vertical coordinates via the TerrainFollowingVerticalDiscretization (TFVD) grid type.

TFVD stores a uniform reference vertical coordinate $r$ and a formulation (e.g. LinearDecay or TwoLevelDecay) that defines the physical altitude

\[z(x, y, r) = r + h(x, y) \, b(r) ,\]

where $h(x, y)$ is the terrain and $b(r)$ is a decay basis satisfying $b(0) = 1$ and $b(z_\text{top}) = 0$.

Public API:

• TerrainFollowingVerticalDiscretization — pass as the z argument to RectilinearGrid / LatitudeLongitudeGrid.
• LinearDecay (Gal-Chen & Somerville 1975).
• TwoLevelDecay (Schär et al. 2002).
• materialize_terrain! — fill the formulation's terrain fields from a h(x, y) function once the grid is built.
• build_terrain_metrics — attach a pressure-gradient stencil.

See docs/src/terrain_following_coordinates.md for the math, discrete operators, well-balancing reference state, and worked examples.

struct LinearDecay{FT, H, SX, SY} <: Breeze.TerrainFollowingDiscretization.AbstractTerrainFormulation

Gal-Chen & Somerville (1975) terrain-following formulation: a single decay basis $b(r) = 1 - r/z_{top}$ that linearly attenuates the terrain from the surface to the model top.

struct SlopeInsideInterpolation

Terrain pressure gradient stencil where the slope is multiplied inside the interpolation of $∂p'/∂r$:

\[\text{correction} = \overline{\overline{s \, \partial_r p'}^x}^z\]

The slope is evaluated at each (Center, Center, Face) stencil point before averaging to (Face, Center, Center).

struct SlopeOutsideInterpolation

Terrain pressure gradient stencil where the slope is multiplied outside the interpolation of $∂p'/∂r$:

\[\text{correction} = s(i,j,k) \, \overline{\overline{\partial_r p'}^x}^z\]

This is the default stencil.

TerrainFollowingVerticalDiscretization(r_faces; formulation=LinearDecay())

Skeleton constructor. r_faces is the reference (computational) r face specification — a range, vector, or function — exactly as for a static z grid. The terrain components inside formulation are filled later by materialize_terrain! once the horizontal grid exists.

struct TerrainMetrics{H, SX, SY, FT, PG}

Pre-computed terrain derivative fields and model top height.

Fields

• topography: 2D field storing $h(x, y)$ at (Center, Center)
• ∂x_h: 2D field storing $\partial h / \partial x$ at (Face, Center)
• ∂y_h: 2D field storing $\partial h / \partial y$ at (Center, Face)
• z_top: Height of the model top (top of the reference coordinate)
• pressure_gradient_stencil: Stencil type for the terrain-corrected horizontal pressure gradient (SlopeOutsideInterpolation or SlopeInsideInterpolation)

struct TwoLevelDecay{ZT, FT, H, SX, SY, B} <: Breeze.TerrainFollowingDiscretization.AbstractTerrainFormulation

Schär et al. (2002) "Smooth LEvel VErtical" (SLEVE) terrain-following formulation. Splits the terrain into a smoothed large-scale component $h_1$ (decay length large_scale_height) and the residual small-scale component $h_2$ (decay length small_scale_height). Each is attenuated with a hyperbolic-sine basis $b_n(r) = \sinh((z_{top}-r)/s_n) / \sinh(z_{top}/s_n)$, so the small-scale features decay quickly while the large-scale envelope is preserved aloft.

Constructed via the kwarg form TwoLevelDecay(; large_scale_height, small_scale_height).

build_terrain_metrics(grid, stencil) -> TerrainMetrics

Build a TerrainMetrics for a materialized TerrainFollowingVerticalDiscretization grid. On such grids the terrain slope used by the dynamics comes from the grid ∂z∂x operator (formulation decay), so this object only carries the pressure_gradient_stencil, z_top, and a representative terrain field.

materialize_terrain!(grid, topography) -> Any

Fill the terrain components of a TerrainFollowingVerticalDiscretization grid in place from topography(x, y). Must be called after the grid is built (the horizontal nodes are needed to evaluate the topography). For TwoLevelDecay, the terrain is split into large- and small-scale parts by horizontal smoothing.

The topography is evaluated at the horizontal cell-centre nodes. Its arguments follow the grid's horizontal coordinates with Flat dimensions dropped, matching every other set! initialiser: topography(x, y) on a RectilinearGrid, topography(λ, φ) on a LatitudeLongitudeGrid, and e.g. topography(x) when y is Flat.