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.

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

Return CondensedPhase with specified parameters converted to FT.

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

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

Arguments

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

struct IdealGas{FT}

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

Fields

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

Examples

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

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

struct MoistureMassFractions{FT}

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

Fields

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

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

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

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

Arguments

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

Keyword arguments

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

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

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

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

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

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

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

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

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

The latent heat varies linearly with temperature:

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

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

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

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

The latent heat varies linearly with temperature:

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

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

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

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

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

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

where:

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

Arguments

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

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

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

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

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

relative_humidity(T, ρ, qᵛ, constants, surface=PlanarLiquidSurface())

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

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

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

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

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

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

Examples

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

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

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

# output
0.010359995391195264

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

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

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

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

# output
0.01128386068542303

saturation_vapor_pressure(T, constants, surface) -> Any

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

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

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

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

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

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

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

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

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

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

which then becomes

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

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

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

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

Arguments

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

vapor_pressure(T, ρ, qᵛ, constants)

Compute the vapor pressure from the ideal gas law:

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

MoistAirBuoyancy(
    grid;
    surface_pressure,
    reference_potential_temperature,
    standard_pressure,
    thermodynamic_constants
) -> MoistAirBuoyancy{RS, AT} where {RS<:(ReferenceState{_A, F} where {_A, F<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B})}), AT<:ThermodynamicConstants}

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

Example

using Breeze, Oceananigans

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

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

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

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

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

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

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

Arguments

• The default dynamics is AnelasticDynamics.

• The default formulation is :LiquidIcePotentialTemperature.

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

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

Example

julia> using Breeze

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

julia> model = AtmosphereModel(grid)
AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=288.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵗ: Centered(order=2)
├── 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.

AtmosphereModelBuoyancy{D, F, T}

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

RadiativeTransferModel(
    grid,
    constants,
    optical_thickness;
    kw...
)

Return a RadiativeTransferModel on grid with thermodynamic constants and using the optical_thickness model for radiative transfer.

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

liquid_ice_potential_temperature(model::AtmosphereModel)

Return an AbstractField representing potential temperature θ for model.

liquid_ice_potential_temperature_density(model::AtmosphereModel)

Return an AbstractField representing potential temperature density for model.

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.

precipitation_rate(model, phase=:liquid)

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

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

Arguments:

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

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

regularize_atmosphere_model_boundary_conditions(boundary_conditions, grid, formulation, thermodynamic_constants)

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

static_energy(model::AtmosphereModel)

Return an AbstractField representing the (specific) static energy for model.

static_energy_density(model::AtmosphereModel)

Return an AbstractField representing static energy density for model.

surface_precipitation_flux(model)

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

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

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

Arguments:

• model: An AtmosphereModel with a microphysics scheme

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

Dynamics

Submodule defining the dynamical formulations for atmosphere models.

Currently supports:

• AnelasticDynamics: filters acoustic waves by assuming density and pressure are small perturbations from a dry, hydrostatic, adiabatic reference state.

Future implementations may include:

• CompressibleDynamics: fully compressible dynamics with prognostic density and pressure.

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

Return AnelasticDynamics representing incompressible fluid dynamics expanded about reference_state.

default_dynamics(grid, constants)

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

default_dynamics(
    grid,
    constants
) -> AnelasticDynamics{R, Nothing} where R<:(ReferenceState{_A, F} where {_A, F<:(Field{Nothing, Nothing, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B})})

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

dynamics_density(dynamics)

Return the density field appropriate to the dynamical formulation.

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

dynamics_density(dynamics::AnelasticDynamics) -> Any

Return the reference density field for AnelasticDynamics.

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

dynamics_pressure(dynamics)

Return the pressure field appropriate to the dynamical formulation.

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

dynamics_pressure(dynamics::AnelasticDynamics) -> Any

Return the reference pressure field for AnelasticDynamics.

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

dynamics_pressure_solver(dynamics, grid)

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

materialize_dynamics(dynamics_stub, grid, boundary_conditions)

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

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

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

materialize_momentum_and_velocities(dynamics, grid, boundary_conditions)

Create momentum and velocity fields for the given dynamics.

mean_pressure(dynamics)

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

mean_pressure(dynamics::AnelasticDynamics) -> Any

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

pressure_anomaly(dynamics)

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

pressure_anomaly(dynamics::AnelasticDynamics) -> Any

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

total_pressure(dynamics)

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

total_pressure(dynamics::AnelasticDynamics) -> Any

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

ThermodynamicFormulations

Submodule defining the thermodynamic formulations for atmosphere models.

A thermodynamic formulation specifies the prognostic thermodynamic variable used to evolve the atmospheric state. Currently supports:

• StaticEnergyFormulation: uses moist static energy density ρe as prognostic variable.
• LiquidIcePotentialTemperatureFormulation: uses liquid-ice potential temperature density ρθ as prognostic variable.

Future implementations may include:

• TotalEnergyFormulation: uses total energy density ρE (for compressible dynamics).

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

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

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

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

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

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

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

additional_thermodynamic_field_names(formulation)

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

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

Collect all prognostic fields into a single NamedTuple.

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

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

compute_thermodynamic_tendency!(model, common_args)

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

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

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

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

diagnose_thermodynamic_state(
    i,
    j,
    k,
    grid,
    formulation::LiquidIcePotentialTemperatureFormulation,
    dynamics,
    q
) -> Breeze.Thermodynamics.LiquidIcePotentialTemperatureState

Build a LiquidIcePotentialTemperatureState at grid point (i, j, k) from the given formulation, dynamics, and pre-computed moisture mass fractions q.

diagnose_thermodynamic_state(
    i,
    j,
    k,
    grid,
    formulation::StaticEnergyFormulation,
    dynamics,
    q
) -> Breeze.Thermodynamics.StaticEnergyState

Build a StaticEnergyState at grid point (i, j, k) from the given formulation, dynamics, and pre-computed moisture mass fractions q.

materialize_formulation(formulation, dynamics, grid, boundary_conditions)

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

Valid symbols:

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

prognostic_thermodynamic_field_names(formulation)

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

thermodynamic_density(formulation)

Return the thermodynamic density field for the given formulation.

thermodynamic_density_name(formulation)

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

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 diagnosing moist instabilities. It is the temperature that a parcel would have if all its moisture were condensed out and the resulting latent heat used to warm the parcel, followed by adiabatic expansion to a reference pressure.

We use a formulation derived by Emanuel (1994),

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

where $T$ is temperature, $pᵈ$ is dry air pressure, $p₀$ is the reference pressure, $ℒˡ$ is the latent heat of vaporization, $qᵛ$ is the vapor specific humidity, $ℋ$ is the relative humidity, and $cᵖᵐ$ is the heat capacity of the moist air mixture.

The formulation follows equation (34) of the paper by Bryan and Fritsch (2002), adapted from the derivation in the work by Durran and Klemp (1982).

Arguments

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

Examples

using Breeze

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

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

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

References

• Bryan, G. H. and Fritsch, J. M. (2002). A benchmark simulation for moist nonhydrostatic numerical models. Monthly Weather Review 130, 2917–2928.
• 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.
• Emanuel, K. A. (1994). Atmospheric Convection (Oxford University Press).

LiquidIcePotentialTemperature(model, flavor=:specific)

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

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

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

or

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

where $θ$ is the potential temperature, $Π = (p/p₀)^{Rᵐ/cᵖᵐ}$ is the Exner function using mixture properties, $ℒˡᵣ$ and $ℒⁱᵣ$ are the reference latent heats of vaporization and sublimation, $qˡ$ and $qⁱ$ are the liquid and ice specific humidities, and $cᵖᵐ$ is the moist air heat capacity.

This is the prognostic thermodynamic variable used in LiquidIcePotentialTemperatureThermodynamics.

Arguments

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

Examples

using Breeze

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

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

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

PotentialTemperature(model, flavor=:specific)

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

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

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

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

Arguments

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

Examples

using Breeze

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

θ = PotentialTemperature(model)
Field(θ)

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

StabilityEquivalentPotentialTemperature(model, flavor=:specific)

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

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

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

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

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

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

Arguments

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

Examples

using Breeze

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

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

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

References

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

StaticEnergy(model, flavor=:specific)

Return a KernelFunctionOperation representing moist static energy $e$.

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

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

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

This is the prognostic thermodynamic variable used in StaticEnergyThermodynamics.

Arguments

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

Examples

using Breeze

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

e = StaticEnergy(model)
Field(e)

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

VirtualPotentialTemperature(model, flavor=:specific)

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

Virtual potential temperature is the temperature that dry air would need to have in order to have the same density as moist air at the same pressure. It accounts for the effect of water vapor on air density:

\[θᵛ = θˡⁱ \left( qᵈ + ε qᵛ \right)\]

where $θˡⁱ$ is liquid-ice potential temperature, $qᵈ$ and $qᵛ$ are the specific humidities of dry air and vapor respectively, and $ε = Rᵛ / Rᵈ$ is the ratio between the vapor and dry air gas constants. $ε ≈ 1.608$ for water vapor and a dry air mixture typical to Earth's atmosphere.

using Breeze

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

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

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

References

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

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

Return a BulkMicrophysics microphysics scheme.

Keyword arguments

• categories: Precipitation categories (e.g., 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

BulkMicrophysics{N, C, B}

Bulk microphysics scheme with cloud formation and precipitation categories.

Fields

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

MixedPhaseEquilibrium(; ...) -> MixedPhaseEquilibrium
MixedPhaseEquilibrium(
    FT;
    freezing_temperature,
    homogeneous_ice_nucleation_temperature
) -> MixedPhaseEquilibrium

Return MixedPhaseEquilibrium representing a temperature-dependent equilibrium between water vapor, possibly supercooled liquid water, and ice.

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

NonEquilibriumCloudFormation(liquid, ice=nothing)

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

The condensation/evaporation rate follows Morrison and Milbrandt (2015), relaxing toward saturation with timescale τ_relax.

Fields

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

References

• Morrison, H., and J. A. Milbrandt (2015). Parameterization of cloud microphysics based on the prediction of bulk ice particle properties. Part I: Scheme description and idealized tests. J. Atmos. Sci., 72, 287–311. https://doi.org/10.1175/JAS-D-14-0065.1

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

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

The options for equilibrium are:

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

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

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

Return a KernelFunctionOperation representing the specified flavor of saturation specific humidity $qᵛ⁺$ which correpsonds to model.microphysics. If model.microphysics is not a saturation adjustment scheme, then a warm phase scheme is assumed which computes the saturation specific humidity over a planar liquid surface.

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 saturated conditions, using the model.specific_moisture. This is equivalent to the :total_moisture flavor under saturated conditions with no condensate; or in other words, if model.specific_moisture happens to be equal to the saturation specific humidity.

• :total_moisture

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

Examples

using Breeze
grid = RectilinearGrid(size=(1, 1, 128), extent=(1e3, 1e3, 1e3))
microphysics = SaturationAdjustment()
model = AtmosphereModel(grid; microphysics)
set!(model, θ=300)
qᵛ⁺ = SaturationSpecificHumidity(model, :prognostic)

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

As SaturationSpecificHumidity it may be wrapped in Field to store the result of its computation. For example, a Field representing the equilibrium saturation specific humidity may be formed via,

qᵛ = SaturationSpecificHumidity(model, :equilibrium) |> Field

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

We also provide a constructor and type alias for the Field itself. For example, to build a Field representing the saturation specific humidity in the case that the total specific moisture is exactly at saturation,

qᵗ = SaturationSpecificHumidityField(model, :total_moisture)

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

Return WarmPhaseEquilibrium representing an equilibrium between water vapor and liquid water.

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

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

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

where $pᵛ$ is the vapor pressure (partial pressure of water vapor) computed from the ideal gas law $pᵛ = ρ qᵛ Rᵛ T$, and $pᵛ⁺$ is the saturation vapor pressure.

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

Examples

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

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

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

ℋ_field = RelativeHumidity(model) |> Field

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

We also provide a convenience constructor for the Field:

ℋ_field = RelativeHumidityField(model)

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

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

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

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

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

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

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

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

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

Fields

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

The horizontal average is computed automatically during update_state!.

Example

using Breeze

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

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

model = AtmosphereModel(grid; forcing)

model.forcing.ρθ

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

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

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

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

Arguments

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

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

Example

using Breeze

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

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

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

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

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

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

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

Keyword Arguments

• direction: The direction of the momentum component (XDirection() or YDirection()). If nothing, the direction is inferred from the field location during boundary condition regularization.
• coefficient: The drag coefficient (default: 1e-3).
• gustiness: Minimum wind speed to prevent singularities when winds are calm (default: 0)

BulkSensibleHeatFluxFunction(; coefficient, gustiness=0, surface_temperature)

Create a bulk sensible heat flux function for computing surface heat fluxes. The flux is computed as:

\[Jᵀ = - ρ₀ Cᵀ |U| (θ - θ₀)\]

where Cᵀ is the transfer coefficient, |U| is the wind speed, θ is the atmospheric potential temperature at the surface, and θ₀ is the surface temperature.

Keyword Arguments

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

BulkVaporFluxFunction(; coefficient, gustiness=0, surface_temperature)

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

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

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

Keyword Arguments

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

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

Create a FluxBoundaryCondition for surface momentum drag.

See BulkDragFunction for details.

Examples

using Breeze

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

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

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

using Oceananigans.Grids: XDirection

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

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

and similarly for YDirection for v.

BulkSensibleHeatFlux(; coefficient, gustiness=0, surface_temperature)

Create a FluxBoundaryCondition for surface sensible heat flux.

See BulkSensibleHeatFluxFunction for details.

Example

using Breeze

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

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

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

BulkVaporFlux(; coefficient, surface_temperature, gustiness=0)

Create a FluxBoundaryCondition for surface moisture flux.

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

See BulkVaporFluxFunction for details.

Example

using Breeze

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

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

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

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

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

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

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

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

The solar zenith angle $θ_z$ satisfies:

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

where:

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

Arguments

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

Returns

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

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

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

equation_of_time(day_of_year) -> Any

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

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

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

References

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

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

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

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

Arguments

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

solar_declination(day_of_year) -> Any

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

Uses the approximation by Spencer (1971):

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

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

References

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

compute_boussinesq_adjustment_temperature(
    𝒰₀::Breeze.Thermodynamics.LiquidIcePotentialTemperatureState{FT},
    constants
) -> Any

Return the temperature $T$ corresponding to thermodynamic equilibrium between the specific humidity and liquid mass fractions of the input thermodynamic state 𝒰₀, wherein the specific humidity is equal to or less than the saturation specific humidity at the given conditions and affiliated with theromdynamic constants constants.

The saturation equilibrium temperature satisfies the nonlinear relation

\[θ = [1 - ℒˡᵣ qˡ / (cᵖᵐ T)] T / Π ,\]

with $ℒˡᵣ$ the latent heat at the reference temperature $Tᵣ$, $cᵖᵐ$ the mixture specific heat, $Π$ the Exner function, $qˡ = \max(0, qᵗ - qᵛ⁺)$ the condensate specific humidity, $qᵗ$ is the total specific humidity, and $qᵛ⁺$ is the saturation specific humidity.

The saturation equilibrium temperature is thus obtained by solving $r(T) = 0$, where

\[r(T) ≡ T - θ Π - ℒˡᵣ qˡ / cᵖᵐ .\]

Solution of $r(T) = 0$ is found via the secant method.

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

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

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

compute_auxiliary_variables!(model)

Compute auxiliary model variables:

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

• thermodynamic variables from the prognostic thermodynamic state,

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

compute_forcings!(model)

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

compute_moisture_fractions(
    i,
    j,
    k,
    grid,
    microphysics::Nothing,
    ρ,
    qᵗ,
    μ
) -> Breeze.Thermodynamics.MoistureMassFractions

Build and return MoistureMassFractions at (i, j, k) for the given grid, microphysics, microphysical_fields, and total moisture mass fraction qᵗ.

Dispatch is provided for ::Nothing microphysics here. Specific microphysics schemes may extend this method to provide tailored behavior.

Note: while ρ and qᵗ are scalars, the microphysical fields μ are NamedTuple of Field. This may be changed in the future.

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

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

maybe_adjust_thermodynamic_state(
    i,
    j,
    k,
    state,
    _::Nothing,
    ρᵣ,
    microphysical_fields,
    qᵗ,
    thermo
) -> Any

Possibly apply saturation adjustment. If a microphysics scheme does not invoke saturation adjustment, just return the state unmodified. In contrast to adjust_thermodynamic_state, this function ingests the entire microphysics formulation and the microphysical_fields. This is needed because some microphysics schemes apply saturation adjustment to a subset of the thermodynamic state (for example, omitting precipitating species).

Grid indices (i, j, k) are provided to allow access to prognostic microphysical fields at the current grid point. The reference density ρᵣ is passed to avoid recomputing it.

microphysical_tendency(
    i,
    j,
    k,
    grid,
    microphysics::Nothing,
    name,
    ρ,
    μ,
    𝒰,
    constants
) -> Any

Return the tendency of the microphysical field name associated with microphysics and thermodynamic constants.

TODO: add the function signature when it is stable

microphysical_velocities(
    microphysics::Nothing,
    microphysical_fields,
    name
)

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

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

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

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

set_thermodynamic_variable!(
    model::AtmosphereModel{<:Any, <:LiquidIcePotentialTemperatureFormulation},
    _::Val{:T},
    value
)

Set the thermodynamic state from in-situ temperature $T$.

The temperature is converted to liquid-ice potential temperature θˡⁱ using the relation between $T$ and θˡⁱ` that accounts for the moisture distribution.

For unsaturated air (no condensate), this simplifies to $θ = T / Π$ where $Π$ is the Exner function.

set_thermodynamic_variable!(
    model::AtmosphereModel{<:Any, <:StaticEnergyFormulation},
    _::Val{:T},
    value
)

Set the thermodynamic state from temperature $T$.

The temperature is converted to static energy $e$ using the relation:

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

settable_specific_microphysical_names(microphysics)

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

specific_to_density_weighted(name::Symbol)

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

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

update_microphysical_fields!(
    microphysical_fields,
    microphysics::Nothing,
    i,
    j,
    k,
    grid,
    density,
    state,
    thermo
)

Update microphysical fields for microphysics_scheme given the thermodynamic state and thermodynamic parameters.

update_radiation!(radiation, model)

Update the radiative fluxes from the current model state.

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

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

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

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

Set variables in an AtmosphereModel.

Keyword Arguments

Variables are set via keyword arguments. Supported variables include:

Prognostic variables (density-weighted):

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

Settable thermodynamic variables:

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

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

• u, v, w: velocity components (sets both velocity and momentum)
• qᵗ: total specific moisture (sets both specific and density-weighted moisture)

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

• qᶜˡ: specific cloud liquid (sets ρqᶜˡ = ρᵣ * qᶜˡ)
• qʳ: specific rain (sets ρqʳ = ρᵣ * qʳ)
• Other prognostic microphysical variables with the ρ prefix removed

Options

• enforce_mass_conservation: If true (default), applies a pressure correction to ensure the velocity field satisfies the anelastic continuity equation.

compute_flux_bc_tendencies!(model::AtmosphereModel)

Apply boundary conditions by adding flux divergences to the right-hand-side.

make_pressure_correction!(
    model::Breeze.AtmosphereModels.AnelasticModel,
    Δt
)

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

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

regularize_atmosphere_model_boundary_conditions(boundary_conditions, grid, surface_pressure, thermodynamic_constants)

Regularize boundary conditions for AtmosphereModel. This function walks through all boundary conditions and calls regularize_atmosphere_boundary_condition on each one, allowing specialized handling for bulk flux boundary conditions and other atmosphere-specific boundary condition types.

RadiativeTransferModel(
    grid,
    constants,
    optical_thickness::RRTMGP.AtmosphericStates.GrayOpticalThicknessOGorman2008;
    surface_temperature,
    coordinate,
    epoch,
    stefan_boltzmann_constant,
    avogadro_number,
    latitude,
    surface_emissivity,
    direct_surface_albedo,
    diffuse_surface_albedo,
    surface_albedo,
    solar_constant
)

Construct a gray atmosphere radiative transfer model for the given grid.

Keyword Arguments

• latitude: Latitude in degrees. If nothing (default), extracted from grid y-coordinate. Can be a scalar (constant for all columns) or a 2D field.
• surface_temperature: Surface temperature in Kelvin (default: 300). Can be a scalar or 2D field.
• surface_emissivity: Surface emissivity, 0-1 (default: 0.98). Scalar.
• direct_surface_albedo: Direct surface albedo, 0-1 (default: 0.1). Can be scalar or 2D field.
• diffuse_surface_albedo: Diffuse surface albedo, 0-1 (default: 0.1). Can be scalar or 2D field.
• solar_constant: Top-of-atmosphere solar flux in W/m² (default: 1361)

update_radiation!(
    rtm::RadiativeTransferModel{<:RRTMGP.AtmosphericStates.GrayOpticalThicknessOGorman2008},
    model
)

Update the radiative fluxes from the current model state.

This function:

1. Updates the RRTMGP atmospheric state from model fields (T, p)
2. Computes the solar zenith angle from the model clock and grid location
3. Solves the longwave and shortwave RTE
4. Copies the fluxes to Oceananigans fields for output

Sign convention: positive flux = upward, negative flux = downward.

copy_fluxes_to_fields!(
    rtm::RadiativeTransferModel{<:RRTMGP.AtmosphericStates.GrayOpticalThicknessOGorman2008},
    grid
)

Copy RRTMGP flux arrays to Oceananigans ZFaceFields.

Applies sign convention:

• positive = upward
• negative = downward.

For the non-scattering shortwave solver, only the direct beam flux is computed.

rrtmgp_context(arch::CPU) -> Any

Create an RRTMGP-compatible ClimaComms context from an Oceananigans architecture.

update_rrtmgp_state!(
    rrtmgp_state::RRTMGP.AtmosphericStates.GrayAtmosphericState,
    model,
    surface_temperature
)

Update the RRTMGP GrayAtmosphericState arrays from model fields.

Grid staggering: layers vs levels

RRTMGP requires atmospheric state at both "layers" (cell centers) and "levels" (cell faces). This matches the finite-volume staggering used in Oceananigans:

                        ┌─────────────────────────────────────────────────┐
    z_lev[Nz+1] ━━━━━━━ │  level Nz+1 (TOA):  p_lev, t_lev, z_lev         │ ← extrapolated
                        └─────────────────────────────────────────────────┘
                        ┌─────────────────────────────────────────────────┐
                        │  layer Nz:  T[Nz], p_lay[Nz] = pᵣ[Nz]           │ ← from model
                        └─────────────────────────────────────────────────┘
    z_lev[Nz]   ━━━━━━━   level Nz:   p_lev, t_lev, z_lev                   ← interpolated
                        ┌─────────────────────────────────────────────────┐
                        │  layer Nz-1                                     │
                        └─────────────────────────────────────────────────┘
                                            ⋮
                        ┌─────────────────────────────────────────────────┐
                        │  layer 2                                        │
                        └─────────────────────────────────────────────────┘
    z_lev[2]    ━━━━━━━   level 2:    p_lev, t_lev, z_lev                   ← interpolated
                        ┌─────────────────────────────────────────────────┐
                        │  layer 1:   T[1], p_lay[1] = pᵣ[1]              │ ← from model
                        └─────────────────────────────────────────────────┘
    z_lev[1]    ━━━━━━━   level 1 (surface, z=0):  p_lev = p₀, t_lev      │ ← from reference state
                        ══════════════════════════════════════════════════
                                        GROUND (t_sfc)

Why the model must provide level values

RRTMGP is a general-purpose radiative transfer solver that operates on columns of atmospheric data. It does not interpolate from layers to levels internally because:

1. Boundary conditions: The surface (level 1) and TOA (level Nz+1) require boundary values that only the atmospheric model knows. For pressure, we use the reference state's surface_pressure at z=0. For the top, we extrapolate using the adiabatic hydrostatic formula.

2. Physics-appropriate interpolation: Different quantities need different interpolation methods. Pressure uses geometric mean (log-linear interpolation) because it varies exponentially with height. Temperature uses arithmetic mean.

3. Model consistency: The pressure profile must be consistent with the atmospheric model's reference state. RRTMGP has no knowledge of the anelastic approximation or the reference potential temperature θ₀.

Physics notes

Temperature: We use the actual temperature field T from the model state. This is the temperature that matters for thermal emission and absorption.

Pressure: In the anelastic approximation, pressure perturbations are negligible compared to the hydrostatic reference pressure. We use reference_state.pressure at cell centers, computed via adiabatic_hydrostatic_pressure(z, p₀, θ₀).

RRTMGP array layout

• Layer arrays (nlay, ncol): values at cell centers, layer 1 at bottom
• Level arrays (nlev, ncol): values at cell faces, level 1 at surface (z=0)
• nlev = nlay + 1

update_solar_zenith_angle!(
    sw_solver,
    coordinate::Tuple,
    grid,
    datetime
)

Update the solar zenith angle in the shortwave solver from the model clock.

Uses the datetime from clock.time and the grid's location (latitude/longitude) to compute the cosine of the solar zenith angle via celestial mechanics.

Does not support anything but single-column grids for now.

OneMomentCloudMicrophysics(FT = Oceananigans.defaults.FloatType;
                           cloud_formation = NonEquilibriumCloudFormation(CloudLiquid(FT), nothing),
                           categories = one_moment_cloud_microphysics_categories(FT),
                           precipitation_boundary_condition = nothing)

Return a OneMomentCloudMicrophysics microphysics scheme for warm-rain and mixed-phase precipitation.

The one-moment scheme uses CloudMicrophysics.jl 1M processes:

• Condensation/evaporation of cloud liquid (relaxation toward saturation)
• Autoconversion of cloud liquid to rain
• Accretion of cloud liquid by rain
• Terminal velocity for rain sedimentation

By default, non-equilibrium cloud formation is used, where cloud liquid is a prognostic variable that evolves via condensation/evaporation tendencies following Morrison and Milbrandt (2015). The prognostic variables are ρqᶜˡ (cloud liquid mass density) and ρqʳ (rain mass density).

For equilibrium (saturation adjustment) cloud formation, pass:

using Breeze.Microphysics
cloud_formation = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())

# output
SaturationAdjustment{WarmPhaseEquilibrium, Float64}(0.001, Inf, WarmPhaseEquilibrium())

Keyword arguments

• precipitation_boundary_condition: Controls whether precipitation passes through the bottom boundary.
  • nothing (default): Rain exits through the bottom (open boundary)
  • ImpenetrableBoundaryCondition(): Rain collects at the bottom (zero terminal velocity at surface)

See the CloudMicrophysics.jl documentation for details.

References

• Morrison, H., and J. A. Milbrandt (2015). Parameterization of cloud microphysics based on the prediction of bulk ice particle properties. Part I: Scheme description and idealized tests. J. Atmos. Sci., 72, 287–311. https://doi.org/10.1175/JAS-D-14-0065.1

ZeroMomentCloudMicrophysics(FT = Oceananigans.defaults.FloatType;
                            cloud_formation = SaturationAdjustment(FT),
                            τ_precip = 1000,
                            qc_0 = 5e-4,
                            S_0 = 0)

Return a ZeroMomentCloudMicrophysics microphysics scheme for warm-rain precipitation.

The zero-moment scheme removes cloud liquid water above a threshold at a specified rate:

• τ_precip: precipitation timescale in seconds (default: 1000 s)

and either

• S_0: supersaturation threshold (default: 0)
• qc_0: cloud liquid water threshold for precipitation (default: 5×10⁻⁴ kg/kg)

For more information see the CloudMicrophysics.jl documentation.

ZeroMomentBulkMicrophysics

Type alias for BulkMicrophysics with CloudMicrophysics 0M precipitation scheme.

The 0M scheme instantly removes precipitable condensate above a threshold. Interface is identical to non-precipitating microphysics except that maybe_adjust_thermodynamic_state calls CloudMicrophysics remove_precipitation first.

precipitation_rate(model, microphysics::ZeroMomentCloudMicrophysics, ::Val{:liquid})

Return a Field representing the liquid precipitation rate (rain rate) in kg/kg/s.

For zero-moment microphysics, this is the rate at which cloud liquid water is removed by precipitation: -dqᵗ/dt from the remove_precipitation function.

surface_precipitation_flux(
    model,
    microphysics::BulkMicrophysics{<:Any, <:Breeze.Microphysics.FourCategories{<:CloudMicrophysics.Parameters.CloudLiquid, <:CloudMicrophysics.Parameters.CloudIce, <:CloudMicrophysics.Parameters.Rain, <:CloudMicrophysics.Parameters.Snow, <:CloudMicrophysics.Parameters.CollisionEff, <:CloudMicrophysics.Parameters.Blk1MVelType, <:CloudMicrophysics.Parameters.AirProperties}}
) -> Field{LX, LY, Nothing, O, G, I, D, T, B, Oceananigans.Fields.FieldStatus{Float64}} where {LX, LY, O, G, I, D, T, B}

Return a 2D Field representing the precipitation flux at the bottom boundary.

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

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

Note: The returned value is positive when rain is falling out of the domain (the terminal velocity wʳ is negative, and we flip the sign).

condensation_rate(qᵛ, qᵛ⁺, qᶜˡ, T, ρ, q, τᶜˡ, constants)

Compute the condensation/evaporation rate for cloud liquid water.

Returns the rate of change of cloud liquid mass fraction (kg/kg/s). Positive values indicate condensation, negative values indicate evaporation.

The rate follows Morrison and Milbrandt (2015); Eq. (A1):

\[\frac{dqᶜˡ}{dt} = \frac{qᵛ - qᵛ⁺}{Γˡ τˡ}\]

Evaporation is limited to the available cloud liquid to prevent negative values.

References

• Morrison, H., and J. A. Milbrandt (2015). Parameterization of cloud microphysics based on the prediction of bulk ice particle properties. Part I: Scheme description and idealized tests. J. Atmos. Sci., 72, 287–311. https://doi.org/10.1175/JAS-D-14-0065.1

deposition_rate(qᵛ, qᵛ⁺ⁱ, qᶜⁱ, T, ρ, q, τᶜⁱ, constants)

Compute the deposition/sublimation rate for cloud ice.

Returns the rate of change of cloud ice mass fraction (kg/kg/s). Positive values indicate deposition, negative values indicate sublimation.

diffusional_growth_factor(aps::AirProperties, T, constants)

Compute the thermodynamic factor $G$ that controls the rate of diffusional growth of cloud droplets and rain drops.

The $G$ factor combines the effects of thermal conductivity and vapor diffusivity on phase change. It appears in the Mason equation for droplet growth:

\[\frac{dm}{dt} = 4π r G 𝒮\]

where $𝒮$ is supersaturation and $r$ is droplet radius.

This is a translation of CloudMicrophysics.Common.G_func_liquid using Breeze's thermodynamics instead of Thermodynamics.jl.

See Eq. (13.28) by Pruppacher & Klett (2010).

References

• Pruppacher, H. R., Klett, J. D. (2010). Microphysics of clouds and precipitation. Springer Netherlands. 2nd Edition

ice_thermodynamic_adjustment_factor(qᵛ⁺ⁱ, T, q, constants)

Compute the thermodynamic adjustment factor $Γⁱ$ for deposition/sublimation.

Same as thermodynamic_adjustment_factor but uses ice latent heat and saturation over ice surface.

rain_evaporation(rain_params, vel, aps, q, qʳ, ρ, T, constants)

Compute the rain evaporation rate (dqʳ/dt, negative for evaporation).

This is a translation of CloudMicrophysics.Microphysics1M.evaporation_sublimation that uses Breeze's internal thermodynamics instead of Thermodynamics.jl.

Arguments

• rain_params: Rain microphysics parameters (pdf, mass, vent)
• vel: Terminal velocity parameters
• aps: Air properties (kinematic viscosity, vapor diffusivity, thermal conductivity)
• q: MoistureMassFractions containing vapor, liquid, and ice mass fractions
• qʳ: Rain specific humidity
• ρ: Air density
• T: Temperature
• constants: Breeze ThermodynamicConstants

Returns

Rate of change of rain specific humidity (negative = evaporation)

thermodynamic_adjustment_factor(qᵛ⁺, T, q, constants)

Compute the thermodynamic adjustment factor $Γˡ$ for condensation/evaporation.

This factor accounts for the temperature dependence of saturation vapor pressure during phase change, following Morrison and Milbrandt (2015); eq. (A1),

\[Γˡ = 1 + \frac{ℒˡ}{cᵖᵐ} \frac{dqᵛ⁺}{dT}\]

where the temperature derivative of saturation specific humidity is:

\[\frac{dqᵛ⁺}{dT} = qᵛ⁺ \left( \frac{ℒˡ}{Rᵛ T²} - \frac{1}{T} \right)\]

References

• Morrison, H., and J. A. Milbrandt (2015). Parameterization of cloud microphysics based on the prediction of bulk ice particle properties. Part I: Scheme description and idealized tests. J. Atmos. Sci., 72, 287–311. https://doi.org/10.1175/JAS-D-14-0065.1