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.

bottom_advective_tracer_flux(i, j, grid, advection, ρ, w, c)

Compute the downward advective mass flux per unit area at the bottom face (k = 1) for tracer c advected by vertical velocity w through density field ρ, using the given advection scheme.

For an explicit scheme this evaluates the same face flux that div_ρUc uses at the bottom boundary. For adaptive implicit vertical advection it returns the instantaneous split-operator flux at the supplied tracer state; the step-integrated implicit outflow depends on the post-solve tracer and must be accumulated separately.

For adaptive implicit vertical advection, the result includes both the CFL-scaled high-order explicit flux and its first-order implicit remainder, evaluated at c.

Returns a positive value for downward (out-of-domain) flux.

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;
    pressure_solver_float_type=Float64
) -> AnelasticDynamics{<:Any, Nothing, DataType}

Return AnelasticDynamics representing incompressible fluid dynamics expanded about reference_state.

Keyword arguments

• pressure_solver_float_type: floating-point type of the anelastic pressure solve. Default: Float64, also for a Float32 model. The Fourier–tridiagonal solve is poorly conditioned for the longest horizontal waves under a tall, strongly stratified column, and can be unstable in Float32 on large domains. When it differs from the grid's type, the solver is built on a copy of the grid in this type, and the pressure is returned in the model's precision.

total_energy_density_name

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

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

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

total_moisture_density_name

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

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

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

AbstractMicrophysicalState{FT}

Abstract supertype for microphysical state structs.

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

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

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

See also microphysical_state, microphysical_tendency.

abstract type AbstractNegativeMoistureCorrection

Abstract supertype for negative moisture correction schemes.

See fix_negative_moisture! for details.

abstract type AbstractNumberConcentrationCategories

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

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

abstract type AbstractSolarPosition

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

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

struct AdiabaticBalancer{T, S}

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

Keyword arguments

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

struct AllSkyOptics <: Breeze.AtmosphereModels.AbstractOptics

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

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=nothing,
    epoch=nothing
) -> 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=Clock(grid),
    thermodynamic_constants=ThermodynamicConstants(eltype(grid)),
    formulation=:LiquidIcePotentialTemperature,
    dynamics=nothing,
    velocities=nothing,
    moisture_density=DefaultValue(),
    tracers=tuple(),
    coriolis=nothing,
    boundary_conditions=NamedTuple(),
    forcing=NamedTuple(),
    advection=DefaultValue(),
    momentum_advection=DefaultValue(),
    scalar_advection=DefaultValue(),
    closure=nothing,
    microphysics=nothing,
    timestepper=nothing,
    timestepper_kwargs=NamedTuple(),
    radiation=nothing,
    particles::Union{Nothing, Oceananigans.TimeSteppers.AbstractLagrangianParticles}=nothing
) -> AtmosphereModel{Dyn, Frm, Arc, Tst, Grd, Clk, Thm, Mom, Moi, Nothing, Tmp, Sol, Vel, Trc, Adv, Nothing, Frc, Nothing, Cnd, Sed, Nothing, Nothing, Nothing, Nothing} where {Dyn, Frm, Arc, Tst, Grd, Clk, Thm, Mom, Moi, Tmp, Sol, Vel, Trc, Adv, Frc, Cnd, Sed}

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

Arguments

• The default dynamics is AnelasticDynamics.

• The default formulation is :LiquidIcePotentialTemperature.

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

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

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

Example

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

References

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

struct AtmosphereModelBuoyancy{D, F, T, M, MF}

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. The microphysics and its fields are carried so that the moisture mass fractions — condensate included — can be recovered at any grid point through grid_moisture_fractions.

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₂=0.78084,
    O₂=0.20946,
    CO₂=0.00042,
    CH₄=1.8e-6,
    N₂O=3.3e-7,
    CO=0.0,
    NO₂=0.0,
    O₃=standard_ozone_profile,
    CFC₁₁=0.0,
    CFC₁₂=0.0,
    CFC₂₂=0.0,
    CCl₄=0.0,
    CF₄=0.0,
    HFC₁₂₅=0.0,
    HFC₁₃₄ₐ=0.0,
    HFC₁₄₃ₐ=0.0,
    HFC₂₃=0.0,
    HFC₃₂=0.0
) -> 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)

CellAdvectionTimescale{F}(
    formulation
) where F -> CellAdvectionTimescale
CellAdvectionTimescale(
    formulation::F
) where F -> CellAdvectionTimescale

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.

The vertical term uses the largest transport speed among the air and every advected sedimenting prognostic (air velocity plus sedimentation velocity), so ThreeDimensionalFormulation() reports the true three-dimensional CFL whether or not the vertical advection is implicit.

struct ClearSkyOptics <: Breeze.AtmosphereModels.AbstractOptics

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

struct CloudScatteringTables{L, I}

Cloud single-scattering tables for EcCKDOptics: one table per condensed phase, each either a Symbol naming a table shipped with the ecRad data (:mie_droplet for spherical liquid droplets, :baum_general_habit_mixture for the Baum general-habit ice mixture) or a path to a netCDF file in the same format.

Fields:

• liquid: Liquid droplet scattering table, a Symbol selector or a file path
• ice: Ice particle scattering table, a Symbol selector or a file path

julia> using Breezejulia> CloudScatteringTables()CloudScatteringTables(liquid=:mie_droplet, ice=:baum_general_habit_mixture)

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(
    FT::DataType=Oceananigans.defaults.FloatType;
    latitude,
    declination=0,
    day_length=86400,
    noon_offset=0
)

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)

EcCKDOptics(
    gas_model=:climate_32x32;
    clouds=nothing,
    lookup_search=nothing
) -> EcCKDOptics{<:Any, Nothing, Nothing}

Construct EcCKDOptics with the reference gas_model (default :climate_32x32) and, for all-sky radiation, CloudScatteringTables as clouds (default nothing: clear sky). lookup_search (default nothing: NumericalRadiation's default) selects how the gas optics tables are bracketed.

struct EcCKDOptics{M, C, S} <: Breeze.AtmosphereModels.AbstractOptics

Full-spectrum radiation with the ecCKD correlated-k gas optics of Hogan & Matricardi (2022), solved column by column by NumericalRadiation.jl. Can be used as the optics argument of RadiativeTransferModel once the extension is loaded with using NumericalRadiation: NumericalRadiation and using NCDatasets (NumericalRadiation exports its own ThermodynamicConstants, so it is loaded qualified next to using Breeze).

The solved column reaches beyond the top of the grid through a OverlyingAtmosphereColumn, so that the downwelling fluxes at the top of the domain include the emission and absorption of the atmosphere above it.

Fields:

• gas_model: Gas optics model: a Symbol or String selecting a reference ecCKD model shipped with the ecRad data (:climate_32x32, :climate_64x64, ...), a (longwave = path, shortwave = path) pair of ecCKD definition files, or a preloaded NumericalRadiation gas optics model
• clouds: CloudScatteringTables for all-sky radiation, or nothing for clear sky
• lookup_search: the NumericalRadiation.AbstractLookupSearch used to bracket the temperature and Planck-source tables of a gas model read from files, or nothing for NumericalRadiation's default. NumericalRadiation.FixedBinarySearch() has no data-dependent loop exit, which Reactant's while-loop raising needs to differentiate through the gas optics

julia> using Breezejulia> EcCKDOptics()EcCKDOptics├── gas_model: :climate_32x32└── clouds: nothingjulia> EcCKDOptics(:climate_64x64, clouds = CloudScatteringTables())EcCKDOptics├── gas_model: :climate_64x64└── clouds: CloudScatteringTables(liquid=:mie_droplet, ice=:baum_general_habit_mixture)

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

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. Solved by RRTMGP.jl's gray solver, which requires using RRTMGP, ClimaComms, NCDatasets.

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.

HydrostaticallyBalancedDensity(; surface_pressure = nothing)

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

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

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

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

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

NothingMicrophysicalState{FT}

A microphysical state with no prognostic variables.

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

OverlyingAtmosphereColumn(
    FT=Oceananigans.defaults.FloatType;
    top=65000.0,
    layers=40,
    blending_height=1000.0,
    temperature=standard_atmosphere_temperature,
    specific_humidity=standard_specific_humidity,
    ozone_mole_fraction=nothing
) -> OverlyingAtmosphereColumn{<:Any, typeof(standard_atmosphere_temperature), typeof(Breeze.AtmosphereModels.standard_specific_humidity), Nothing, Nothing, Nothing}

Construct a OverlyingAtmosphereColumn of float type FT reaching top (default 65 km) with layers layers (default 40), a blending_height of 1 km for the temperature anchor, and the standard_atmosphere_temperature and standard_specific_humidity profiles.

struct OverlyingAtmosphereColumn{FT, T, Q, O, B, V}

The atmosphere above the top of the grid, as seen by a column radiation solver.

Radiation is solved on the grid's Nz layers plus layers extension layers stacked from the grid's top face z_top to top, geometrically stretched from a first layer as thick as the grid's top layer: with ratio r between successive layers, Δz₁ (rᴺ - 1) / (r - 1) = top - z_top (see overlying_atmosphere_faces). No layers are added when the grid already reaches top.

The extension temperature follows temperature(z), anchored to the temperature T_top on the grid's top face so that the profile is continuous there and relaxes to temperature(z) over blending_height,

\[Tₑ(z) = temperature(z) + (T_top - temperature(z_top)) exp(-(z - z_top) / blending_height) ,\]

with blending_height = 0 disabling the anchor. Pressure is hydrostatic above the top face, the specific humidity follows specific_humidity(z), and the ozone mole fraction follows ozone_mole_fraction(z), or the O₃ of the BackgroundAtmosphere when ozone_mole_fraction is nothing. Each profile is a number (constant with height) or a function of height. Extension layers are clear.

A radiation model materializes the extension on its grid: the materialized OverlyingAtmosphereColumn also holds the extension's layers and the profiles sampled at their centers, which are nothing until then.

Fields:

• top: Height of the top of the extended column [m]
• layers: Number of extension layers between the grid top and top
• blending_height: Height over which the temperature anchor at the grid top decays [m]; 0 disables it
• temperature: Temperature profile [K]
• specific_humidity: Specific humidity profile [kg kg⁻¹]
• ozone_mole_fraction: Ozone mole fraction profile [mol mol⁻¹], or nothing for BackgroundAtmosphere.O₃
• bottom: Height of the bottom of the extension, the grid's top face [m]
• profile_temperature_at_bottom: The temperature profile at bottom [K], from which the anchor measures the grid top's departure
• Δzᵃᵃᶜ: Thickness of the extension layers, bottom-up [m]
• zᵃᵃᶜ: Height of the extension layers' centers [m]
• T: The temperature profile at the layer centers, before anchoring [K]
• qᵛ: The specific_humidity profile at the layer centers [kg kg⁻¹]
• O₃: The ozone mole fraction profile at the layer centers [mol mol⁻¹]

julia> using Breezejulia> OverlyingAtmosphereColumn()OverlyingAtmosphereColumn{Float64}├── top: 65000.0 m├── layers: 40├── blending_height: 1000.0 m├── temperature: standard_atmosphere_temperature (generic function with 1 method)├── specific_humidity: standard_specific_humidity (generic function with 1 method)└── ozone_mole_fraction: nothingjulia> OverlyingAtmosphereColumn(Float32, top = 80e3, layers = 60, ozone_mole_fraction = 1e-6)OverlyingAtmosphereColumn{Float32}├── top: 80000.0 m├── layers: 60├── blending_height: 1000.0 m├── temperature: standard_atmosphere_temperature (generic function with 1 method)├── specific_humidity: standard_specific_humidity (generic function with 1 method)└── ozone_mole_fraction: 1.0e-6

RadiativeTransferModel(
    grid::Oceananigans.Grids.AbstractGrid,
    optics,
    constants=nothing;
    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
• EcCKDOptics() - Full-spectrum clear- or all-sky radiation using ecCKD gas optics via NumericalRadiation.jl (requires using NumericalRadiation: NumericalRadiation and using NCDatasets)

The RRTMGP optics require using RRTMGP, ClimaComms, NCDatasets. NumericalRadiation exports its own ThermodynamicConstants, so load it qualified (using NumericalRadiation: NumericalRadiation) next to using Breeze.

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> using NumericalRadiation: NumericalRadiationjulia> 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)julia> RadiativeTransferModel(grid, EcCKDOptics(), 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)├── liquid_effective_radius: ConstantRadiusParticles{Float64}(1.0e-5)├── ice_effective_radius: ConstantRadiusParticles{Float64}(3.0e-5)├── diffuse_surface_albedo: ConstantField(0.1)├── optics: EcCKDOptics with 32 longwave and 32 shortwave g-points, clear sky└── overlying_atmosphere: 40 layers from 10000.0 m to 65000.0 m

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 SlowDynamics{D}

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

When computing momentum tendencies with a SlowDynamics-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]

MicrophysicalProcessRate(
    model,
    name::Symbol
) -> KernelFunctionOperation{<:Any, <:Any, <:Any, <:Any, T, K, D} where {T, K<:Breeze.Microphysics.PredictedParticleProperties.P3ProcessRateKernelFunction, D<:Tuple}

Return a KernelFunctionOperation at (Center, Center, Center) that evaluates the microphysical process rate called name of model.microphysics (for example :autoconversion, :accretion, :condensation, or :rain_evaporation).

Rates are specific: per unit mass of moist air, in kg kg⁻¹ s⁻¹ for mass transfers and kg⁻¹ s⁻¹ for number transfers. Multiply by the density (for example total_density(model.dynamics)) for a rate per unit volume, whose vertical Integral is a column-integrated rate in kg m⁻² s⁻¹.

Each scheme documents which rates it provides and their sign conventions. Schemes that compute their rates jointly with sink limiters evaluate the limited rates, that is, exactly the rates that drive the model tendencies, at the current model state. The operation is lazy: wrap it in Field to compute and store it, slice it with Field(op; indices), or reduce it with Integral and Average.

Supported schemes: PredictedParticlePropertiesMicrophysics, for which every rate field of P3ProcessRates is available.

Example

using Breezegrid = RectilinearGrid(size=(4, 4, 8), extent=(1e3, 1e3, 1e3))model = AtmosphereModel(grid; microphysics=P3Microphysics())set!(model, θ=290, qᵛ=0.01, qᶜˡ=1e-3)autoconversion = MicrophysicalProcessRate(model, :autoconversion)# outputKernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: P3ProcessRateKernelFunction(:autoconversion,)└── arguments: ()

advecting_momentum(model)

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

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

aerosol_field_names(microphysics) -> Tuple{}

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

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

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

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

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

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.

bottom_precipitation_flux(
    model
) -> Field{LX, LY, LZ, O, G, I, D, T, B, Oceananigans.Fields.FieldStatus{Float64}} where {LX, LY, LZ, O, G, I, D, T, B}

Return a 2D Field representing the flux of precipitating moisture at the bottom boundary, in kg/m²/s (positive = downward flux out of domain).

The default sums the bottom-face advective flux of every sedimenting condensate of the model (see materialize_sedimentation), each evaluated with the same advection scheme that transports that tracer, so the diagnostic matches the boundary flux the tendency operator applies. For a scheme whose cloud condensate also sediments (the non-equilibrium 1M and 2M schemes, and P3) this includes cloud droplet and cloud ice fallout, not only rain and snow. For adaptive implicit vertical advection, this is the split-operator flux evaluated at the current state rather than a time-integrated mass-loss diagnostic. Schemes that move precipitation by internal means (such as DCMIP2016KM) extend bottom_precipitation_flux(model, microphysics) instead.

Building the returned Field assembles a kernel operation, so construct it once and compute! it repeatedly rather than calling this inside a callback.

Arguments:

• model: An AtmosphereModel with a microphysics scheme

The generic method is implemented in Breeze.Advection, which owns the flux kernel and is loaded after this module.

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)

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

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

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

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

compute_pressure_correction!(model, Δt)

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

condensate_liquid_fraction(microphysics, ::Val)

Return the liquid fraction of the condensate mass name: 1 for liquid, 0 for ice, a number in between for a mass of mixed composition, or nothing (the default) for anything that is not a sedimenting condensate mass.

The fraction must match how moisture_fractions bins the mass into the liquid and ice mass fractions, since those fix the latent heat and heat capacity the thermodynamic state attributes to it. Every name in condensate_field_names with a sedimentation_velocity must declare it; materialize_sedimentation checks this at model construction, since a falling mass without one would leave its latent heat behind. The fraction is the enthalpy the mass carries, not what it falls with: P3's liquid on ice ρqʷⁱ falls at the ice speed yet is liquid, because no fusion enthalpy has been released for it and moisture_fractions counts it in the liquid mass fraction.

The default nothing covers two cases, neither needing a method: condensate that is not a prognostic mass (saturation-adjusted cloud), and a tracer that is not a condensate mass (a number moment such as ρnᶜˡ, or a non-additive property such as P3's ρqᶠ and ρbᶠ). A sedimenting condensate mass with no liquid fraction is an error.

The acoustic substep loop of the compressible core projects the tendency of every name in condensate_field_names onto its pressure with the sensitivity of the mass's phase, so there every such name must declare its fraction, whether or not it sediments.

A mixed-phase particle is normally two condensate masses sharing a fall speed, as P3 carries ice ρqⁱ and the liquid on it ρqʷⁱ. One mass of mixed composition may instead return its liquid fraction, which liquid_fraction_content resolves exactly, the content being linear in composition.

continuity_density(dynamics)

Return the continuity density: the density whose continuity equation the model solves, $∇·(ρᵣ u) = 0$ on the anelastic core and $∂ₜρᵈ + ∇·(ρᵈ u) = 0$ on the compressible core. Momentum and the thermodynamic variable are stored per unit of it (ρu = ρᵈ u, ρθ = ρᵈ θ), so it is also the divisor for diagnosing velocity (u = ρu/ρᵈ) and potential temperature (θ = ρθ/ρᵈ).

• 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 quantity, total_density: it recovers the moisture mass fractions (qˣ = ρˣ/ρ) and weights scalar and water advection. The water densities (ρqᵛ, ρqˡ, …) are stored as partial densities (mass per volume), not per unit of the continuity density. On the anelastic core the two densities coincide (total_density === continuity_density).

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 cᵖᵐ T = ĉᵖ (ρ Rᵐ T / pˢᵗ)^κ θ + ℒˡ₀ qˡ + ℒⁱ₀ qⁱ) extend this function to return an iterative solver.

dynamics_pressure(dynamics)

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

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

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

grid_moisture_fractions(
    i,
    j,
    k,
    grid,
    microphysics,
    ρ,
    qᵛ,
    μ_fields
)

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) -> Int64

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

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

Returns 0 by default.

initial_aerosol_number_density(microphysics, ρ)

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

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

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

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

is_density_tendency_forcing(_) -> Bool

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

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

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

liquid_ice_potential_temperature(model)

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

liquid_ice_potential_temperature_density(model)

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

make_pressure_correction!(model, Δt)

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

materialize_atmosphere_model_boundary_conditions(boundary_conditions, grid, formulation,
                                                dynamics, microphysics, thermodynamic_constants,
                                                closure = nothing)

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

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

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

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

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

The turbulence closure sets the gustiness that bulk fluxes constructed without one take: the subgrid turbulent kinetic energy under a TKEBasedTurbulenceClosure, a constant otherwise.

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

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

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

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

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

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

This is called internally by RadiativeTransferModel constructors.

materialize_surface_property(x, grid [, solar_position])

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

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

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

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

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

Arguments

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

Returns

An AbstractMicrophysicalState subtype containing the local specific microphysical variables.

See also microphysical_tendency, AbstractMicrophysicalState.

microphysical_tendencies(
    microphysics,
    names::Tuple,
    ρ,
    ℳ,
    𝒰,
    constants
)

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

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

See also microphysical_tendency, compute_microphysical_tendencies!.

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

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

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

Arguments

• microphysics: The microphysics scheme
• name: Variable name as Val(:name) (e.g., Val(:ρqᶜˡ))
• ρ: Total density, total_density, which weights the density-weighted prognostics and their tendencies (scalar). Because the prognostic condensate content is $ρ q$, it also sets the precipitation size distributions. On the anelastic core this is the reference density $ρᵣ$.
• ℳ: Microphysical state (e.g., WarmPhaseOneMomentState)
• 𝒰: Thermodynamic state
• constants: Thermodynamic constants

Saturation must not use ρ. The saturation specific humidity $qᵛ⁺ = pᵛ⁺ / (ρᵗʰ Rᵛ T)$ and the supersaturation use the thermodynamic state density density(𝒰, constants), which thermodynamic_state_density also returns: $ρᵗʰ = pᵣ / (Rᵐ T)$ on the anelastic core and the total density on the compressible core.

Returns

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

See also microphysical_state, AbstractMicrophysicalState.

microphysical_velocities(
    microphysics,
    microphysical_fields,
    name
) -> Union{Nothing, NamedTuple{(:u, :v, :w), <:Tuple{Oceananigans.Fields.ZeroField{Int64, 3}, Oceananigans.Fields.ZeroField{Int64, 3}, Any}}}

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

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

The generic implementation calls sedimentation_velocity and uses the result as the vertical velocity component.

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_flux(
    model;
    condensates=keys(model.sedimentation)
) -> Union{KernelFunctionOperation{<:Any, <:Any, <:Any, <:Any, T, typeof(Breeze.AtmosphereModels.precipitation_flux_kernel), D} where {T, D<:Tuple{Any, Any}}, KernelFunctionOperation{Center, Center, Face, <:Any, T, typeof(Breeze.AtmosphereModels.precipitation_flux_kernel), D} where {T, D<:Tuple{Any, Any, Any}}}

Return a KernelFunctionOperation at (Center, Center, Face) for the downward mass flux of sedimenting condensate through every $z$-face, in kg m⁻² s⁻¹ (positive downward). Multiply by 3600 to obtain a precipitation rate in mm h⁻¹ of liquid water equivalent.

The flux is the sedimentation part of each condensate's advective mass flux, $-ρ (F(wᵗ + wˢ) - F(wᵗ)) / A_z$, formed with the same advection scheme and face reconstruction that transports the condensate (sedimentation_mass_fluxes plus the adaptive-implicit remainder implicit_sedimentation_mass_fluxes), so it matches the sedimentation the model applies. Here $F(w)$ is the volume-weighted face flux of the condensate's specific humidity at velocity $w$, $wᵗ$ the transport velocity and $wˢ$ the condensate's (negative) sedimentation velocity. At the bottom face, where an impenetrable boundary makes $wᵗ = 0$, it reduces to bottom_precipitation_flux.

The keyword argument condensates selects which sedimenting condensates contribute, by prognostic name (see keys(model.sedimentation)); by default all of them do, which for a scheme with sedimenting cloud condensate includes cloud droplet fallout. Use, for example, condensates = (:ρqʳ,) for the rain flux alone. Sedimentation velocities are those diagnosed during the most recent update_state!.

Example

using Breezegrid = RectilinearGrid(size=(4, 4, 8), extent=(1e3, 1e3, 1e3))model = AtmosphereModel(grid; microphysics=P3Microphysics())set!(model, θ=290, qᵛ=0.01, qᶜˡ=1e-3)P = precipitation_flux(model, condensates=(:ρqʳ,))# outputKernelFunctionOperation at (Center, Center, Face)├── grid: 4×4×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: precipitation_flux_kernel (generic function with 1 method)└── arguments: ("Tuple", "Field", "Field")

precipitation_rate(model, phase::Symbol=:liquid)

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

The precipitation rate is the rate at which moisture is removed from the atmosphere by precipitation processes.

Arguments:

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

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

pressure_anomaly(dynamics)

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

prognostic_field_names(::Nothing) -> Tuple{}

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

prognostic_moisture_phase_fractions(
    i,
    j,
    k,
    grid,
    microphysics,
    microphysical_fields,
    q
) -> Tuple{Any, Any, Any}

Return the fractions (vapor, liquid, ice) of a change of the prognostic moisture density (moisture_prognostic_name) that each phase takes at fixed temperature, in a cell whose moisture mass fractions are q.

The default is all vapor: a scheme with prognostic condensate carries the vapor as its moisture density. A scheme whose moisture density includes condensate that it diagnoses (saturation adjustment) extends this with the partition its equilibrium gives a change.

The acoustic substep loop of the compressible core projects the moisture tendency onto its pressure with the sensitivities of these phases.

resolve_column_batch_rows(::Nothing, Ny)

Number of whole j-rows of columns in each radiation batch when the domain's Ny rows are split into at most column_batches batches. nothing means one unbatched solve over all Ny rows.

Batches are whole rows so that each is a contiguous run of columns, and cld(Ny, column_batches) rows each so that every batch has the same width. When those rows do not divide Ny, the domain needs fewer than column_batches batches.

saturation_adjusted_equilibrium(microphysics)

Return the phase equilibrium to which the scheme adjusts the partition of its prognostic moisture density (moisture_prognostic_name) in saturated air, or nothing if it adjusts none.

The default is nothing: a scheme with prognostic condensate holds the partition it carries. Saturation adjustment returns its equilibrium. In a cell where such a scheme holds cloud condensate (where prognostic_moisture_phase_fractions gives the vapor no share), the pressure that the compressible core re-diagnoses is that of the adjusted partition, in which warming evaporates condensate, and the acoustic substep loop linearizes that pressure.

sedimentation_composition_increment(
    dynamics,
    q,
    phase
) -> Breeze.Thermodynamics.MoistureMassFractions

Return the change of the moisture mass fractions per unit mass of condensate phase that sediments out of a cell whose mass fractions are q, as a MoistureMassFractions difference vector (the dry-air change is implied by Σq = 1). Continuity advances the continuity density (continuity_density) without a sedimentation source on every core, so the increment depends on what the total density does:

• fixed (the default; AnelasticDynamics): the total density the mass fractions are referenced to does not respond, so the dry mass fraction qᵈ = 1 − qᵛ − qˡ − qⁱ makes up the departed mass and the increment is q̂ˣ − q̂ᵈ, with q̂ˣ the composition with all mass in x;
• diagnosed and falling with the condensate (CompressibleDynamics: prognostic dry density ρᵈ, ρ = ρᵈ + Σρˣ): every mass fraction renormalizes and the increment is q̂ˣ − q.

condensate_content differentiates the thermodynamic variable along this increment, ∂φ/∂qˣ = ∇_q φ · Δq: the content per unit falling mass at which the cell the condensate leaves keeps its temperature.

sedimentation_velocity(
    microphysics,
    microphysical_fields,
    ::Val
)

Return the sedimentation velocity field (vertical component, [m/s], negative = downward) for the prognostic tracer name, or nothing (the default) if the tracer does not sediment.

Microphysics schemes extend this function for each sedimenting tracer, dispatching on ::Val{name}.

sedimentation_velocity_field(
    grid
) -> Field{Center, Center, Face, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}

Build a ZFaceField suitable for storing a sedimentation velocity: bottom = nothing ensures a kernel-set bottom-face value is preserved during fill_halo_regions!, while the default impenetrable top holds w = 0 so nothing falls in through the model top.

set_to_mean!(ref::ExnerReferenceState, model)

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

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

set_to_mean!(
    ref::ReferenceState,
    model;
    rescale_densities=false
)

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

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

specific_humidity(model)

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(microphysics, qᵗ, μ, ρ)

Convert total specific moisture $qᵗ$ to the scheme-dependent specific moisture $qᵛᵉ$ given the density-weighted microphysical variables μ and the total air density ρ,

\[qᵛᵉ = qᵗ - Σ ρqᶜ / ρ\]

where the sum runs over the independent condensate mass densities named by condensate_field_names. Diagnostic cloud fractions, number concentrations, and dependent moments (such as the P3 rime mass $ρqᶠ$, already contained in $ρqⁱ$) are therefore excluded. Both the values of μ and ρ may be scalars or fields. For Nothing microphysics there is no condensate, so $qᵛᵉ = qᵗ$.

This is used both by set!, to convert a total moisture input, and by parcel models, which carry total moisture $qᵗ$ as their prognostic variable.

Since condensate_field_names is the complement of the prognostic moisture — the same partition total_condensate_density uses — this recovers vapor for non-equilibrium schemes and equilibrium moisture $qᵉ$ for saturation-adjustment schemes, with no per-scheme method needed.

using Breezeusing Breeze.Microphysics: DCMIP2016KesslerMicrophysicsmicrophysics = DCMIP2016KesslerMicrophysics()μ = (ρqᶜˡ=0.0012, ρqʳ=0.0024)specific_prognostic_moisture(microphysics, 0.02, μ, 1.2)# output0.017

specific_prognostic_moisture(model)

Return the prognostic specific moisture field for model.

This is $qᵛ$ for non-equilibrium schemes or $qᵉ$ for saturation adjustment schemes. To convert total specific moisture instead, pass microphysics, qᵗ, the density-weighted microphysical variables, and the total air density.

specific_thermodynamic_field(formulation)

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

standard_atmosphere_temperature(z)

The temperature (K) of the U.S. Standard Atmosphere 1976 (NOAA, NASA & USAF (1976)) at height z (m): piecewise linear through the seven layers between 0 and 86 km, with the troposphere cooling at 6.5 K km⁻¹ from 288.15 K at the surface to 216.65 K at 11 km, an isothermal lower stratosphere, warming to 270.65 K at the stratopause (47–51 km), and cooling again through the mesosphere to 186.946 K at 84.852 km. Held constant above that, and z is taken as the geopotential height.

julia> using Breezejulia> standard_atmosphere_temperature(0)288.15julia> standard_atmosphere_temperature(11e3)216.65julia> standard_atmosphere_temperature(50e3)270.65

standard_ozone_profile(z)

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.

thermodynamic_density(formulation)

Return the thermodynamic density field for the given formulation — the prognostic thermodynamic variable in continuity-density-weighted ("flux") form (ρθ, ρs, ρE).

The weighting density is the dynamics' continuity density (see continuity_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 θ = ρθ / continuity_density(dynamics).

thermodynamic_density_name(formulation)

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

total_density(
    i,
    j,
    k,
    dry_density,
    microphysics,
    moisture_density,
    microphysical_fields
)

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

total_density(dynamics)

Return the model's total air density ρ = ρᵈ + Σρˣ: the density of the model's mass, which converts the stored partial densities to mass fractions ($qˣ = ρqˣ / ρ$) and weights scalar advection, hydrostatics and radiation. CompressibleDynamics diagnoses it from the prognostic dry density and the partial densities. Defaults to continuity_density, so on AnelasticDynamics it is the reference density $ρᵣ(z)$: the anelastic approximation replaces the total density by $ρᵣ$, consistent with the reference pressure $pᵣ$.

Thermodynamic relations of the local state, such as the saturation specific humidity, use thermodynamic_state_density instead.

total_pressure(dynamics)

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

transport_velocities(model)

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.

validate_schedule(schedule, arch)

Return schedule if radiation can be updated on it on architecture arch, and throw an ArgumentError otherwise. Every schedule is valid on native architectures; BreezeReactantExt accepts only IterationInterval on ReactantState. AtmosphereModel validates the schedule of its radiation when it is constructed, before anything is traced.

write_sedimentation_velocity!(w_field, i, j, k, bc, 𝕎)

Store the fall-speed magnitude 𝕎 diagnosed at cell center (i, j, k) in the sedimentation velocity field w_field at the cell's bottom face: microphysics libraries return positive magnitudes while Breeze stores signed vertical velocities (negative = downward), and sedimentation is always downward, so the donor cell for face k is cell k itself. At k = 1 the precipitation boundary condition bc is applied (see bottom_sedimentation_velocity); the top face (k = Nz + 1) lies outside the :xyz launch region and is held at zero by the impenetrable top boundary condition.

DewpointTemperature(
    model;
    solver=SecantSolver(eltype(model.grid); reltol = 0.0001, abstol = 0, maxiter = 10)
) -> KernelFunctionOperation{<:Any, <:Any, <:Any, <:Any, T, K, D} where {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 StaticStridesOffsetArray(::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 StaticStridesOffsetArray(::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 invariant under phase changes at fixed pressure and total water, and accounts for the latent heat associated with liquid water and ice:

\[θˡⁱ = \frac{cᵖᵐ T - ℒˡ₀ qˡ - ℒⁱ₀ qⁱ}{ĉᵖ Π}\]

where $ĉᵖ = qᵈ cᵖᵈ + qᵗ cᵖᵛ$ is the vapor-equivalent heat capacity, $ℒˡ₀$ and $ℒⁱ₀$ are the latent heats extrapolated to absolute zero, $qˡ$, $qⁱ$ are the liquid and ice mass fractions, and $Π$ is the Exner function with the exponent of model.formulation (see Breeze.PotentialTemperatureFormulations.LiquidIcePotentialTemperatureFormulation). This is the variable that LiquidIcePotentialTemperatureFormulation evolves.

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 StaticStridesOffsetArray(::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 StaticStridesOffsetArray(::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,
    flavor_symbol=:prognostic
) -> KernelFunctionOperation{<:Any, <:Any, <:Any, <:Any, T, K, D} where {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 StaticStridesOffsetArray(::Array{Float64, 3}, 0:2, 0:2, -2:11) with eltype Float64 with indices 0:2×0:2×-2:11    └── max=326.162, min=325.851, mean=326.006

References

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

StaticEnergy(model, flavor=:specific)

Return a KernelFunctionOperation representing moist static energy $s$.

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

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

where $cᵖᵐ$ is the moist air heat capacity, $T$ is temperature, $g$ is gravitational acceleration, $z$ is height, and $ℒˡ₀ qˡ + ℒⁱ₀ qⁱ$ uses latent-heat offsets extrapolated to zero kelvin.

This is the prognostic thermodynamic variable used in StaticEnergyThermodynamics.

Arguments

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

Examples

using Breezegrid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1e3))model = AtmosphereModel(grid)set!(model, θ=300)s = StaticEnergy(model)Field(s)# output1×1×8 Field{Center, Center, Center} on RectilinearGrid on CPU├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo├── boundary conditions: FieldBoundaryConditions│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing├── operand: KernelFunctionOperation at (Center, Center, Center)├── status: time=0.0└── data: 3×3×14 StaticStridesOffsetArray(::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 StaticStridesOffsetArray(::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=(0, 0), m=4)

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=(0, 0), m=4)

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=0.001,
    gustiness=DefaultGustiness(),
    surface_temperature=nothing,
    filtered_velocities=nothing
) -> Breeze.BoundaryConditions.BulkDragFunction{Nothing, Nothing, Float64, Breeze.BoundaryConditions.DefaultGustiness, Nothing, Nothing, Nothing}

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

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

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

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

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

Monin–Obukhov consistency

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

Default surface temperature

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

Keyword Arguments

• direction: The direction of the momentum component (XDirection(), YDirection(), or ZDirection()). If nothing, the direction is inferred from the field location during boundary condition regularization.
• coefficient: The drag coefficient (default: 1e-3). Can be a constant or a PolynomialCoefficient for wind and stability-dependent transfer coefficients.
• gustiness: The gustiness Uᵍ added to the resolved wind, in the flux and in a PolynomialCoefficient alike: a number, a TKEGustiness, or the default DefaultGustiness, which the model replaces with a TKEGustiness under a TKEBasedTurbulenceClosure and with 0 otherwise.
• surface_temperature: Surface temperature, used to compute ρˢ and required when using PolynomialCoefficient with stability correction. Can be a Field, Function, or Number. A function takes the non-Flat coordinates of the wall followed by the time, as for Oceananigans boundary conditions: (x, y, t) on the bottom and top, (y, z, t) on the west and east, (x, z, t) on the south and north. (default: nothing)
• filtered_velocities: A FilteredSurfaceVelocities for temporally filtered wind speed, near-surface velocity, and θᵥ in the bulk formula. If nothing (default), instantaneous fields are used.

BulkSensibleHeatFluxFunction(;
    coefficient,
    gustiness=DefaultGustiness(),
    surface_temperature,
    filtered_velocities=nothing
)

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

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

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

• For LiquidIcePotentialTemperatureFormulation: $Δϕ = θ - θˢ$ (potential temperature flux).
• For StaticEnergyFormulation: $Δϕ = s - sˢ$ (static energy flux).

The wall value uses the atmospheric vapor, liquid, and ice fractions at the sampling point, the wall temperature $Tˢ$, and the wall pressure and height. This holds composition fixed when computing sensible heat exchange. When filtering is enabled, the complete difference $Δϕ$ is sampled at the filter height and filtered in time.

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

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

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

Keyword Arguments

• coefficient: The sensible heat transfer coefficient.
• gustiness: The gustiness $Uᵍ$ added to the resolved wind, $|U| = \sqrt{u² + v² + Uᵍ²}$, in the flux and in a PolynomialCoefficient alike: a number, a TKEGustiness, or the default DefaultGustiness, which the model replaces with a TKEGustiness under a TKEBasedTurbulenceClosure and with 0 otherwise.
• surface_temperature: The wall temperature. Can be a Field, a Function, or a Number. Functions are evaluated at the wall at every time step with the non-Flat coordinates of the wall followed by the time, as for Oceananigans boundary conditions: (x, y, t) on the bottom and top, (y, z, t) on the west and east, (x, z, t) on the south and north, and for example (x, t) on the bottom of a grid that is Flat in y.
• filtered_velocities: Either nothing (default) or FilteredSurfaceVelocities. Note that when filtered_velocities is not nothing, then automatically there is filtering in the scalar fields via FilteredSurfaceScalar with the same parameters (e.g., height, timescale) as filtered_velocities. Filtering is supported on the bottom boundary only.

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

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

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

where $Cᵛ$ is the transfer coefficient, $|U|$ is the wind speed tangential to the wall, $qᵛ$ is the near-wall specific humidity, and $qˢ$ is the specific humidity of the air in contact with the wall. Over the wet fraction $β$ of the wall (the moisture_availability) that is the saturation specific humidity $qᵛ⁺$ at the wall temperature $Tˢ$ times the wall relative humidity $ℋˢ$ (unity for a wet wall); over the dry fraction it is the humidity of the air itself, so that $qᵛ - qˢ = β (qᵛ - ℋˢ qᵛ⁺)$ and the flux is $β$ times the flux over a wet wall. When filtering is enabled, the complete difference $qᵛ - ℋˢ qᵛ⁺$ is sampled at the filter height and filtered in time, so the wall humidity and the near-wall humidity are seen at the same times, as for BulkSensibleHeatFluxFunction.

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

Keyword Arguments

• coefficient: The vapor transfer coefficient.
• gustiness: The gustiness $Uᵍ$ added to the resolved wind, $|U| = \sqrt{u² + v² + Uᵍ²}$, in the flux and in a PolynomialCoefficient alike: a number, a TKEGustiness, or the default DefaultGustiness, which the model replaces with a TKEGustiness under a TKEBasedTurbulenceClosure and with 0 otherwise.
• surface_temperature: The wall temperature. Can be a Field, a Function, or a Number. Used to compute the saturation specific humidity at the wall. Functions take the non-Flat coordinates of the wall followed by the time, as for Oceananigans boundary conditions.
• surface_relative_humidity: The relative humidity of the air in contact with the wall, between 0 and 1 (default: 1, a saturated wall). Can be a Field, a Function, or a Number.
• moisture_availability: The fraction $β ∈ [0, 1]$ of the wall that is wet. nothing (default) takes the value carried by a PolynomialCoefficient coefficient, whose stability correction uses the same surface humidity, and 1 (a wet wall, an ocean) for a constant coefficient. A value that disagrees with a PolynomialCoefficient is an error. The phase of the surface water follows the coefficient in the same way, and is liquid for a constant coefficient.
• filtered_velocities: Either nothing (default) or FilteredSurfaceVelocities. Note that when filtered_velocities is not nothing, then automatically there is filtering in the scalar fields via FilteredSurfaceScalar with the same parameters (e.g., height, timescale) as filtered_velocities. Filtering is supported on the bottom boundary only.

DefaultGustiness() -> Breeze.BoundaryConditions.DefaultGustiness

The placeholder gustiness of a bulk-flux boundary condition constructed without one. When the AtmosphereModel materializes its boundary conditions it is replaced by materialize_gustiness: with TKEGustiness when the closure is a TKEBasedTurbulenceClosure, and with a constant gustiness of zero otherwise.

using Breeze.BoundaryConditions: DefaultGustiness, materialize_gustinessmaterialize_gustiness(DefaultGustiness(), nothing)# output0

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 vapor-equivalent heat capacity $ĉᵖ$ and the Exner function $Π$.

The relationship is:

\[Jᶿ = (ρᵈ / ρ) 𝒬ᵀ / (ĉᵖ Π)\]

where $𝒬ᵀ$ is the energy flux, $Jᶿ$ is the flux of $ρθ = ρᵈ θ$, and $ρᵈ / ρ$ (unity for anelastic dynamics) converts a change of enthalpy per unit total mass into one carried by the dry air. At fixed moisture and pressure $ĉᵖ Π θˡⁱ = cᵖᵐ T - ℒˡ₀ qˡ - ℒⁱ₀ qⁱ$ gives $ĉᵖ Π δθ = cᵖᵐ δT$, so an enthalpy input $𝒬ᵀ = ρ cᵖᵐ \overline{w'T'}$ reaches $θ$ as $Jᶿ = 𝒬ᵀ / (ĉᵖ Π)$. The condensate term cancels, so the conversion holds saturated as well as dry. This is the same conversion the $ρE$ interior forcing and the radiative flux divergence apply. In the same spirit, BulkSensibleHeatFlux refers its wall temperature to $θ$, though at the wall pressure rather than the pressure of the adjacent cell.

The heat capacity and Exner function use moisture fractions from the microphysics scheme and the Exner exponent of the formulation (exner_exponent).

FilteredSurfaceScalar(
    grid;
    filter_timescale,
    height=nothing
)

A two-dimensional field storing a temporally filtered near-surface scalar for use in bulk flux boundary conditions. The bulk scalar fluxes store the complete near-wall difference of their scalar ($Δθ$, $Δs$, or $Δqᵛ$), so the wall value is filtered together with the atmospheric value.

The filter update is the same exponential form as FilteredSurfaceVelocities.

Keyword Arguments

• filter_timescale (required): Filter time scale τ in seconds; must be positive and finite.
• height: Reference height (m) for scalar evaluation. If nothing, the first grid cell center value is used.

FilteredSurfaceVelocities(
    grid;
    filter_timescale,
    height=nothing
)

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

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

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

where τ is the filter_timescale. There is no default: τ sets how much of the near-surface variability the wall model sees, so it must be chosen for the flow at hand. For a wall model driven by the instantaneous near-surface state, omit filtered_velocities from the bulk boundary condition instead.

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

Keyword Arguments

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

FittedStabilityFunction(
    scalar_roughness_length;
    richardson_number_mapping=RichardsonNumberMapping(typeof(scalar_roughness_length)),
    stability_function=StabilityFunction(typeof(scalar_roughness_length))
) -> FittedStabilityFunction{<:Any, RM, SP} where {RM<:Breeze.BoundaryConditions.RichardsonNumberMapping, SP<:Breeze.BoundaryConditions.StabilityFunction}

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

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

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

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

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

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

Arguments

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

Keyword Arguments

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

References

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

PolynomialCoefficient(
    FT=Oceananigans.defaults.FloatType;
    polynomial=nothing,
    roughness_length=0.00015,
    minimum_wind_speed=0.1,
    stability_function=FittedStabilityFunction(FT(roughness_length / 7.3)),
    surface=PlanarLiquidSurface(),
    moisture_availability=1,
    transfer_type=nothing
)

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

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

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

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

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

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

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

Keyword Arguments

• polynomial: Tuple (a₀, a₁, a₂) for the polynomial. If nothing, the polynomial is automatically selected by the boundary condition constructor.
• roughness_length: Surface roughness ℓʳ in meters (default: 1.5e-4, typical for ocean)
• minimum_wind_speed: Minimum wind speed to avoid singularity in a₂/U term (default: 0.1 m/s)
• stability_function: Stability correction strategy. Default is FittedStabilityFunction using Li et al. (2010) $Riᴮ → ζ$ mapping with Hogström (1996) / Beljaars & Holtslag (1991) MOST stability functions. The scalar roughness length defaults to roughness_length / 7.3 (typical ocean value). Use nothing to disable stability correction.
• surface: The phase of the surface water, which selects the saturation specific humidity at the surface in the stability correction: PlanarLiquidSurface() (default), PlanarIceSurface() or PlanarMixedPhaseSurface(liquid_fraction).
• moisture_availability: The fraction $β ∈ [0, 1]$ of the surface that is saturated (default: 1, an ocean). The surface specific humidity entering the stability correction is $qˢ = β qᵛ⁺(Tˢ) + (1 - β) qᵛ$, with $qᵛ$ the specific humidity of the air in the first cell, so that $β = 0$ describes a dry surface whose virtual potential temperature carries no moisture contribution of its own. See wall_virtual_potential_temperature.

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

Examples

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

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

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

References

• Beljaars, A. C. M., & Holtslag, A. A. M. (1991). Flux parameterization over land surfaces for atmospheric models. Journal of Applied Meteorology, 30, 327-341.
• Hogström, U. L. F. (1996). Review of some basic characteristics of the atmospheric surface layer. Boundary-Layer Meteorology, 78, 215-246.
• Large, W., & Yeager, S. G. (2009). The global climatology of an interannually varying air–sea flux data set. Climate dynamics, 33(2), 341-364.
• Li, Y., Gao, Z., Lenschow, D. H., & Chen, F. (2010). An improved approach for parameterizing surface-layer turbulent transfer coefficients in numerical models. Boundary-Layer Meteorology, 137, 153-165.

RichardsonNumberMapping(
    FT=Oceananigans.defaults.FloatType;
    stable_unstable_transition=0,
    strongly_stable_transition=0.2,
    aᵘ₁₁=0.045,
    bᵘ₁₁=0.003,
    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.5,
    bʷ₁₁=-0.0539,
    bʷ₁₂=1.54,
    bʷ₂₁=-0.669,
    bʷ₂₂=-3.282,
    aˢ₁₁=0.7529,
    aˢ₂₁=14.94,
    bˢ₁₁=0.1569,
    bˢ₂₁=-0.3091,
    bˢ₂₂=-1.303
) -> Breeze.BoundaryConditions.RichardsonNumberMapping

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

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

Three regimes:

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

References

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

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

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.

TKEGustiness(
    FT=Oceananigans.defaults.FloatType;
    coefficient=sqrt(2 / 3),
    minimum=0.5
) -> TKEGustiness

A gustiness $Uᵍ$ estimated from the subgrid turbulent kinetic energy $e$ of a TKEBasedTurbulenceClosure in the cell next to the wall,

\[Uᵍ = \max(c \sqrt{e}, Uᵍ_{min}),\]

where $c$ is the coefficient and $Uᵍ_{min}$ the minimum. The default $c = \sqrt{2/3}$ is the isotropic estimate: for isotropic turbulence $e = (u'^2 + v'^2 + w'^2) / 2$ and the variance of each velocity component is $(2/3) e$, so $c \sqrt{e}$ is the root-mean-square of one horizontal component of the subgrid wind. The default minimum of 0.5 m/s is the floor Large and Yeager (2009) place on the wind speed in their bulk formulas.

The gustiness enters the bulk fluxes through the effective wind $Ũ = \sqrt{u² + v² + Uᵍ²}$, which both the flux and its transfer coefficient see. It requires the prognostic TKE density ρe of a TKEBasedTurbulenceClosure, and is the gustiness a bulk-flux boundary condition takes by default under one (see DefaultGustiness).

Keyword Arguments

• coefficient: The coefficient $c$ in $Uᵍ = c \sqrt{e}$ (default: $\sqrt{2/3}$).
• minimum: The floor $Uᵍ_{min}$ on the gustiness, in m/s (default: 0.5).

Example

using BreezeTKEGustiness(minimum=0.2)# outputTKEGustiness{Float64}(coefficient=0.816497, minimum=0.2)

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 vapor-equivalent heat capacity $ĉᵖ$ and the Exner function $Π$.

The relationship is:

\[𝒬ᵀ = Jᶿ ĉᵖ Π\]

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

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

Create a FluxBoundaryCondition for wall momentum drag, on any of the six boundaries.

See BulkDragFunction for details.

Examples

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

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

using Oceananigans.Grids: XDirectionu_drag = BulkDrag(direction=XDirection(), coefficient=1e-3)ρu_bcs = FieldBoundaryConditions(bottom=u_drag)# outputOceananigans.FieldBoundaryConditions, with boundary conditions├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)├── south: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)├── north: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)├── bottom: FluxBoundaryCondition: BulkDragFunction(direction=XDirection(), coefficient=0.001, gustiness=DefaultGustiness())├── top: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)

and similarly for YDirection for v. The same condition may be placed on the walls of a closed box; the direction is inferred from the momentum component it is attached to:

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

BulkSensibleHeatFlux(; coefficient, gustiness=DefaultGustiness(), surface_temperature)

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

The bulk formula computes

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

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

See BulkSensibleHeatFluxFunction for details.

Example

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

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

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

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

See BulkVaporFluxFunction for details.

Example

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

EnergyFluxBoundaryCondition(
    flux,
    κ̂=VaporEquivalentExponent()
)

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 vapor-equivalent heat capacity $ĉᵖ$ and the Exner function $Π$, and by $ρ / ρᵈ$, to obtain the potential temperature flux: $Jᶿ = (ρᵈ / ρ) 𝒬ᵀ / (ĉᵖ Π)$. The model constructor replaces κ̂ with the Exner exponent of its formulation.

ThetaFluxBoundaryCondition(
    flux,
    κ̂=VaporEquivalentExponent()
)

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 vapor-equivalent heat capacity $ĉᵖ$ and the Exner function $Π$, with Exner exponent κ̂, to obtain the energy flux: $𝒬ᵀ = (ρ / ρᵈ) Jᶿ ĉᵖ Π$.

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

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

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

cos_solar_zenith_angle(
    datetime::Dates.AbstractDateTime,
    longitude,
    latitude
)

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

The solar zenith angle $θ_z$ satisfies:

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

where:

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

Arguments

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

Returns

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

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

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

equation_of_time(day_of_year)

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

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

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

References

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

hour_angle(datetime::Dates.AbstractDateTime, longitude)

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)

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

Uses the approximation by Spencer (1971):

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

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

References

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

CompressibleEquations

Module implementing fully compressible dynamics for atmosphere models.

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

The fully compressible Euler equations in conservation form are:

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

Pressure is computed from the ideal gas law:

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

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

abstract type AbstractRamp

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

Most shapes have compact support, meaning they are exactly zero below the layer base $z_{\rm sponge\_top} - \text{depth}$; GaussianRamp deliberately does not, and instead carries a weak tail all the way down.

Choose a shape on absorption behaviour alone. UpperSponge evaluates the ramp once per z-face at construction and stores the resulting 1D profile, so the substep loop never sees the shape, however expensive it is to evaluate.

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, TD, MP, TAV, GT, TS, WC, DC, TWC, IS}

Storage and parameters for the split-explicit acoustic substepper (scheme described in the module header). θᴸ=ρθᴸ/ρᴸ and Cᴸ 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).
• minimum_horizontal_spacing: the grid's smallest horizontal spacing, which the adaptive count measures the acoustic Courant number against; computed once, when the substepper is built.
• forward_weight: CN off-centering ω (0.5 = centered; default 0.65).
• damping, substep_distribution: divergence-damping strategy; substep allocation across WS-RK3 stages.
• linearization_potential_temperature (θᴸ), linearization_pressure_coefficient (Cᴸ, the slope ∂p/∂(ρθ) of the equation of state at the stage-entry state): per-stage caches, see refresh_linearization_basic_state!.
• density_perturbation (ρ′), density_potential_temperature_perturbation ((ρθ)′), momentum_perturbation ((ρu/v/w)′ as .u/.v/.w): perturbation prognostics advanced in the loop.
• density_predictor, density_potential_temperature_predictor: explicit predictors before the vertical solve.
• previous_density_potential_temperature_perturbation: prior-substep (ρθ)′, for Klemp 2018 damping.
• time_averaged_velocities: acoustic-mean velocity for non-acoustic scalar transport (moisture/tracers/ chemistry/TKE); the slow ρθ tendency uses the current RK predictor velocity instead, not this cache.
• slow_vertical_momentum_tendency (Gˢρw, z-faces): advection+Coriolis+closure+forcing (PGF/buoyancy excluded — those are in the fast operator).
• vertical_solver_source_term (z-faces): explicit RHS of the (ρw)′ tridiagonal system.
• vertical_solver: factorized BatchedTridiagonalSolver of the implicit (ρw)′ update, factored once per stage (factor_vertical_solver!) so that each substep solves by substitution.
• vertical_velocity_cache, density_cache: the stage-entry predictor velocity and the dry density that weights it, frozen for implicit_substep! (see cache_advecting_state!); nothing without adaptive-implicit advection.
• time_averaged_vertical_velocity_cache: the acoustic-mean transport velocity the moisture and tracer tendencies were built with, frozen for scalar_substep! and for the thermodynamic tendency's condensate sedimentation term (see cache_transport_velocity!).

AcousticSubstepper(
    grid,
    split_explicit::SplitExplicitTimeDiscretization;
    prognostic_momentum=nothing,
    substep_floattype=eltype(grid),
    cache_advecting_state=false
)

Construct an AcousticSubstepper with auxiliary BCs for perturbations and transport velocities. Bounded normal faces default to nothing; the loop explicitly zeros lateral momentum perturbations. Prescribed flow stays in the stage-entry state, with prognostic BCs reapplied after stage recovery. Inheriting momentum BCs would impose full-state targets on perturbations (issue #716) and velocities. The prognostic_momentum kwarg is retained for backwards compatibility but no longer consulted.

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

Fully compressible dynamics with prognostic density and diagnostic pressure.

Fields

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

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

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

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

CompressibleDynamics(
    time_discretization::TD=ExplicitTimeStepping();
    standard_pressure=100000.0,
    base_pressure=101325.0,
    reference_potential_temperature=nothing,
    reference_temperature=nothing,
    reference_vapor_mass_fraction=nothing,
    slope_stencil=SlopeOutsideInterpolation(),
    terrain_metrics=nothing,
    reference_state=:auto,
    temperature_tolerance=nothing,
    temperature_maxiter=nothing,
    surface_pressure=nothing
) where TD -> CompressibleDynamics{<:Any, Nothing, Nothing, Nothing, Float64, Breeze.CompressibleEquations.AutoReference, SlopeOutsideInterpolation, Nothing, Nothing}

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

Positional Arguments

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

Keyword Arguments

• standard_pressure: Reference pressure for potential temperature (default: 10⁵ Pa)

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

• reference_potential_temperature: Potential temperature for building a fixed hydrostatically-balanced reference state used in base-state subtraction. Can be a constant θ₀ or a function θ(z). Default: nothing, which uses the automatic θᵣ = 288 K profile when reference_state = :auto. When provided, this profile replaces the automatic profile when building the ExnerReferenceState.

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

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

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

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

  Deep near-isentropic reference profiles

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

struct ConstantSubstepSize <: AcousticSubstepDistribution

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

struct CubicRamp <: AbstractRamp

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

Has zero derivative at both the layer base and the lid, so absorbs upgoing waves smoothly without the reflective kink of LinearRamp. Functionally equivalent to Sin2Ramp. The default, and the shape to prefer when the dynamics just below the sponge must be provably untouched; see GaussianRamp for the trade against far-field cleanliness.

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

Gaussian sponge ramp, $\exp[-(z - H)^2 / \text{depth}^2]$, reaching 1 at the lid and $e^{-1}$ one depth below it. Unlike CubicRamp, Sin2Ramp and LinearRamp it is not clamped to zero at the layer base: a weak tail continues all the way down.

That tail is the point: on a long mountain-wave integration a compactly-supported ramp leaves a residual in the far field, worst downstream of the ridge, that a tailed one removes at half the depth. Neither deepening nor strengthening the clamped ramp substitutes — both come out worse downstream, and with a slightly larger core wave, so what they add is reflection rather than absorption. A clamped ramp's relative gradient (1/γ)|∂z γ| is unbounded at its base, and deepening it only moves that base closer to the wave; a Gaussian has no base to reflect off.

What the tail costs is weak damping below the nominal layer. For a steadily forced wave that is not exp(-rate × ramp × t): the response equilibrates against advection, so the amplitude cost scales with rate × ramp(z) / (U k), which on the Schär case is 10% at z = H/2 and 0.07% at H/4. Measured rms w over 0.5 < z < 10 km differs from the clamped ramp by 0.4%. Prefer a clamped ramp when the dynamics just below the sponge must be provably untouched, and this one when a quiet far field matters more.

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. Provided for parity with WRF (damp_opt = 3) / MPAS-Atmosphere.

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, M, P}

Implicit upper Rayleigh sponge for the substepper. Damps the vertical momentum 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 damping acts on the total momentum, which the substep loop carries as $ρw = ρw^L + (ρw)′$, so it is applied in two matching pieces:

• the stage-entry part $-\text{rate} × \text{ramp}(z) × ρw^L$ is a known constant across the substeps and enters the slow tendency $G^s_{ρw}$. This half targets the Cartesian $ρw$ on every grid, height-coordinate and terrain-following alike, because a horizontally uniform wind over sloping coordinate surfaces has $ρw = 0$ but $ρ\tilde{w} = -\text{slope} · ρu$, so damping $ρ\tilde{w}^L$ would drive the mean flow to follow the coordinate surfaces;
• the acoustic perturbation part is CN-weighted inside the column tridiag (paralleling the implicit divergence damping): $δτᵐ⁺ × \text{rate} × \text{ramp}(z)$ on the LHS diagonal, $δτˢ⁻ × \text{rate} × \text{ramp}(z)$ on the explicit-half RHS. This half acts on whichever vertical-momentum perturbation the substep evolves, $(ρw)′$ in height coordinates and $(ρ\tilde{w})′$ in terrain-following coordinates. Either way it needs no Cartesian correction, since the mean-flow offset lives entirely in the stage-entry state.

Damping $(ρw)′$ alone would not absorb anything: the perturbation measures the change over an RK stage, so a quasi-steady gravity wave passes through untouched. This split 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(), GaussianRamp(), LinearRamp(). Custom shapes are supported by subtyping AbstractRamp and defining (::MyRamp)(z, sponge_top, depth).

• mask: optional extra damping region, a function mask(x, y, z) of the face node coordinates (λ, φ, z on a LatitudeLongitudeGrid) returning a value in $[0, 1]$. The damping rate becomes $γ = \text{rate} × \max(\text{ramp}(z), \text{mask}(x, y, z))$, so the same implicit Rayleigh treatment covers, for example, the lateral relaxation zone of a nested limited-area model. Default nothing (lid layer only). Because the damping is implicit, it removes the rate × Δt ≲ 2.5 stability limit that an explicit Relaxation forcing on ρw imposes on the outer time step.

Without a mask, 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. That is also what makes the damping rate $γ(z) = \text{rate} × \text{ramp}(z)$ a 1D profile over z-faces: UpperSponge(; ...) is a lightweight skeleton with damping_profile = nothing, and materialize_sponge evaluates ramp once per face at AcousticSubstepper construction and stores the result. The substep kernels then read γ with a single load instead of re-evaluating the shape twice per face per acoustic substep.

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ᴸ,
    advection=nothing;
    rounding=RoundToNearest()
)

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. The recovery of the full state from the stage-entry state and the accumulated perturbations is rounded with the round-off treatment rounding.

compute_transport_density_tendency!(
    substepper::AcousticSubstepper,
    model,
    velocities
)

Compute the dry air's share, ρᵈ / ρ, of the divergence of the mass flux ρ 𝐮 that the moisture tendencies were built with, where velocities are their transport velocities: the dry density tendency that the dilution part of the moisture content ramps with (Section 3b). Call this whenever the moisture tendencies are built, with the velocities they were built with.

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_potential_temperature = θᴸ = ρθᴸ/ρᴸ derived from model.dynamics.dry_density + ρθ
• linearization_pressure_coefficient = Cᴸ, the slope ∂p/∂(ρθ) of the equation of state at that state (see refresh_linearization_basic_state!)

prepare_acoustic_cache!(
    substepper::AcousticSubstepper,
    model
)

Stage-start cache preparation. Refreshes the cached linearization quantities (θᴸ, Cᴸ) 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{<:Any, Nothing, Nothing}

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

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

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

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

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

SubsidenceForcing(
    wˢ;
    advection=UpwindBiased(order = 1)
) -> SubsidenceForcing{<:Any, Nothing, S, Nothing} where S<:(UpwindBiased{N, FT, Oceananigans.Utils.ExplicitTimeDiscretization} where {N, FT})

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.

Arguments

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

Keyword arguments

• advection: the discretization of $w^s ∂_z \overline{ϕ}$. The default UpwindBiased(order = 1) damps $2Δz$ modes and does not undershoot at sharp gradients. Centered() cannot see a $2Δz$ mode, so it persists; higher-order upwind schemes damp it but may over- or undershoot.

With an upwind scheme every cell, including the two boundary cells, is updated from its own face fluxes. Air leaving the column through a boundary face carries the boundary cell's value. Air entering it carries the value extrapolated linearly from the two cells inside, so a linear profile crosses the boundary undistorted; for mass fractions such as moisture that value is bounded below by zero, so subsidence never brings negative moisture into the column.

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└── advection: UpwindBiased(order=1)

TemperatureForcing(
    temperature_tendency
) -> TemperatureForcing{<:Any, Nothing, Nothing, Nothing, Nothing, Nothing}

Forcing by a prescribed temperature tendency $∂ₜT$ (K s⁻¹), such as the large-scale heating of a single-column case. Supply it under the specific energy key E, as in forcing = (; E = TemperatureForcing(∂ₜT)). ∂ₜT can be anything accepted as a forcing: a Field, a function, a Returns, an Oceananigans.Forcing. At kernel time it is converted to an energy tendency with the local mixture heat capacity,

\[F_E = cᵖᵐ ∂ₜT ,\]

where $cᵖᵐ$ is evaluated from the moisture mass fractions the tendency kernels see (vapor, liquid and ice from the microphysics scheme). Breeze then multiplies by the density and routes the result onto the thermodynamic variable: LiquidIcePotentialTemperatureFormulation divides an energy tendency by $cᵖᵐ Π$, so $θ$ changes at $∂ₜT / Π$, the rate at which $T = Π θ + (ℒˡᵣ qˡ + ℒⁱᵣ qⁱ) / cᵖᵐ$ changes at fixed moisture and pressure by $∂ₜT$. StaticEnergyFormulation adds it to $ρs$ unconverted, and $s = cᵖᵐ T + g z - ℒˡᵣ qˡ - ℒⁱᵣ qⁱ$ changes by $cᵖᵐ ∂ₜT$.

Converting with the dry heat capacity $cᵖᵈ$ instead, as in E = cᵖᵈ * ∂ₜT, heats humid air too little by the factor $cᵖᵈ / cᵖᵐ$: about 1.7% at 20 g kg⁻¹ of vapor.

using Breezegrid = RectilinearGrid(size=4, z=(0, 1000), topology=(Flat, Flat, Bounded))cooling(z, t) = -2 / 86400 # radiative cooling, K s⁻¹forcing = TemperatureForcing(cooling)model = AtmosphereModel(grid; forcing=(; E=forcing))forcing# outputTemperatureForcing with temperature_tendency: cooling (generic function with 1 method)

geostrophic_forcings(uᵍ, vᵍ)

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.

PrescribedDensity{D}(density) where D -> PrescribedDensity
PrescribedDensity(density::D) where D -> PrescribedDensity

Wrapper indicating that density is fixed (not prognostic).

struct PrescribedDynamics{Div, D, P, SP, FT}

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

PrescribedDynamics(
    density;
    pressure=nothing,
    surface_pressure=nothing,
    base_pressure=101325,
    standard_pressure=100000.0,
    divergence_correction=false
) -> PrescribedDynamics{<:Any, <:Any, Nothing, Nothing}

Construct PrescribedDynamics from a density field or PrescribedDensity. If pressure=nothing, hydrostatic pressure is computed during materialization. base_pressure is the pressure datum at z = 0. On a raised domain, the default bottom-face pressure extends the lowest prescribed density down to the datum; pass surface_pressure to prescribe a different hydrostatic anchor explicitly.

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

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

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

Example

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

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

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

Return a BulkMicrophysics microphysics scheme.

Keyword arguments

• categories: Microphysical categories (e.g., cloud liquid, cloud ice, rain, snow) or nothing for non-precipitating
• cloud_formation: Cloud formation scheme (default: SaturationAdjustment)
• precipitation_boundary_condition: Bottom boundary condition for precipitation sedimentation.
  • nothing (default): Precipitation passes through the bottom
  • ImpenetrableBoundaryCondition(): Precipitation collects at the bottom
• negative_moisture_correction: Correction scheme for negative moisture produced by advection.
  • nothing (default): No correction
    • VerticalBorrowing(): Vertical redistribution of the moisture prognostic only
  • SpeciesBorrowing(): Same-level species borrowing only
  • SpeciesBorrowing(vertical_borrowing=VerticalBorrowing()): Species borrowing with vertical redistribution

struct BulkMicrophysics{N, C, B, NMC}

Bulk microphysics scheme with cloud formation and precipitation categories.

Fields

• cloud_formation: Cloud formation scheme (saturation adjustment or non-equilibrium)
• categories: Precipitation categories (e.g., rain, snow) or nothing
• precipitation_boundary_condition: Bottom boundary condition for precipitation sedimentation.
  • nothing (default): Precipitation passes through the bottom (open boundary)
  • ImpenetrableBoundaryCondition(): Precipitation collects at the bottom (zero terminal velocity at surface)
• negative_moisture_correction: Correction scheme for negative moisture produced by advection.
  • nothing (default): No correction
    • VerticalBorrowing(): Vertical redistribution of the moisture prognostic only
  • SpeciesBorrowing(): Same-level species borrowing only
  • SpeciesBorrowing(vertical_borrowing=VerticalBorrowing()): Species borrowing with vertical redistribution

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

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

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

DCMIP2016KesslerMicrophysics(
    FT=Oceananigans.defaults.FloatType;
    dcmip_temperature_scale=237.3,
    terminal_velocity_coefficient=36.34,
    density_scale=0.001,
    terminal_velocity_exponent=0.1364,
    autoconversion_rate=0.001,
    autoconversion_threshold=0.001,
    accretion_rate=2.2,
    accretion_exponent=0.875,
    evaporation_ventilation_coefficient_1=1.6,
    evaporation_ventilation_coefficient_2=124.9,
    evaporation_ventilation_exponent_1=0.2046,
    evaporation_ventilation_exponent_2=0.525,
    diffusivity_coefficient=2.55e8,
    thermal_conductivity_coefficient=540000.0,
    substep_cfl=0.8,
    maximum_sedimentation_substeps=0
) -> Breeze.Microphysics.DCMIP2016KesslerMicrophysics