Thermodynamics

Thermodynamics

Moist thermodynamic functions for atmospheric modeling, including air pressure, latent heats, saturation vapor pressures, and saturation specific humidities.

Parameter Sets

Many functions defined in this module rely on ClimaParams.jl. ClimaParams.jl defines several functions (e.g., many planet parameters). For example, to compute the mole-mass ratio of dry air and water vapor:

import ClimaParams as CP
import Thermodynamics.Parameters as TP
FT = Float64
param_set = TP.ThermodynamicsParameters(FT)
Rv_over_Rd = TP.Rv_over_Rd(param_set)

Because these parameters are widely used throughout this module, param_set is an argument for many Thermodynamics functions.

Saturation adjustment

One entry point is saturation_adjustment. It accepts:

• method: A root-solving method type from RootSolvers.jl (e.g., RS.NewtonsMethod, RS.SecantMethod).
• param_set: The thermodynamics parameter set.
• formulation: The thermodynamic formulation (e.g., ρe(), ph()).
• Thermodynamic state variables appropriate for the formulation.

Supported methods in RootSolvers.jl:

• NewtonsMethod: Newton method with analytic gradients (recommended; analytic derivatives exist for all formulations).
• NewtonsMethodAD: Newton method with automatic differentiation.
• SecantMethod: Secant method (derivative-free).
• BrentsMethod: Brent's method (hybrid root-finding).

Thermodynamic functions

Many thermodynamic functions (e.g., air_temperature, air_pressure) are dispatched over the independent variables used to compute them. The available independent variable types (subtypes of IndepVars) are:

• ρe: Density, internal energy, and specific humidities
• pe: Pressure, internal energy, and specific humidities
• ph: Pressure, enthalpy, and specific humidities
• pρ: Pressure, density, and specific humidities
• pθ_li: Pressure, liquid-ice potential temperature, and specific humidities
• ρθ_li: Density, liquid-ice potential temperature, and specific humidities

Parameters

Submodule defining thermodynamic parameter types and accessor functions.

AbstractThermodynamicsParameters{FT}

Abstract supertype for thermodynamic parameter sets.

Subtypes:

• ThermodynamicsParameters: default concrete parameter set.

Subtypes must implement accessor functions for all thermodynamic constants (e.g., T_0, R_d, cp_d, …). See ThermodynamicsParameters for the full list of required fields.

ThermodynamicsParameters{FT} <: AbstractThermodynamicsParameters{FT}

Concrete parameter set for Thermodynamics.jl.

Fields

• T_0: reference temperature [K].
• T_triple: triple-point temperature of water [K].
• T_freeze: freezing temperature of water [K].
• T_icenuc: temperature of homogeneous ice nucleation [K].
• T_min: lower bound of the saturation-adjustment root search [K]; bounds the search, not the returned temperature.
• T_max: maximum allowable temperature for saturation adjustment [K].
• T_init_min: no longer used; retained for backward compatibility [K]. Saturation adjustment starts from the unsaturated temperature directly, floored only by a numerics bound far below any physical temperature.
• T_surf_ref: reference surface temperature [K].
• T_min_ref: minimum reference temperature (for temperature profiles) [K].
• entropy_reference_temperature: reference temperature for entropy [K].
• MSLP: mean sea-level pressure [Pa].
• p_ref_theta: reference pressure for potential temperature [Pa].
• press_triple: triple-point pressure of water [Pa].
• R_d: specific gas constant of dry air [J/(kg·K)].
• R_v: specific gas constant of water vapor [J/(kg·K)].
• cp_d: isobaric specific heat capacity of dry air [J/(kg·K)].
• cp_v: isobaric specific heat capacity of water vapor [J/(kg·K)].
• cp_l: isobaric specific heat capacity of liquid water [J/(kg·K)].
• cp_i: isobaric specific heat capacity of ice [J/(kg·K)].
• LH_v0: latent heat of vaporization at T_0 [J/kg].
• LH_s0: latent heat of sublimation at T_0 [J/kg].
• entropy_dry_air: reference specific entropy of dry air [J/(kg·K)].
• entropy_water_vapor: reference specific entropy of water vapor [J/(kg·K)].
• grav: gravitational acceleration [m/s²].
• pow_icenuc: exponent for power-law ice nucleation fraction [-].
• q_min: minimum condensate threshold for has_condensate [kg/kg].

Examples

import ClimaParams as CP
import Thermodynamics.Parameters as TP

FT = Float64;
param_set = TP.ThermodynamicsParameters(FT)

# Alternatively:

toml_dict = CP.create_toml_dict(FT)
param_set = TP.ThermodynamicsParameters(toml_dict)

Ratio of gas constants R_v / R_d.

Latent heat of fusion at reference temperature T_0.

Internal energy of vaporization at reference temperature T_0.

Internal energy of fusion at reference temperature T_0. Equal to the latent heat of fusion because for incompressible condensed phases the p Δv term is negligible (p Δv ≪ L_f).

Ratio of dry air gas constant to isobaric specific heat R_d / cp_d.

Isochoric specific heat capacity of dry air.

Isochoric specific heat capacity of water vapor.

Isochoric specific heat capacity of liquid water.

Isochoric specific heat capacity of ice.

Phase

A condensed phase, to dispatch over (e.g., in saturation_vapor_pressure).

Subtypes:

• Liquid: liquid water phase.
• Ice: ice phase.

Liquid <: Phase

A liquid phase, to dispatch over (e.g., in saturation_vapor_pressure).

Ice <: Phase

An ice phase, to dispatch over (e.g., in saturation_vapor_pressure).

IndepVars

Abstract type to dispatch over the independent variables used to compute a thermodynamic quantity (e.g., in air_temperature).

Subtypes:

• ρe: density, internal energy, and specific humidities.
• pe: pressure, internal energy, and specific humidities.
• ph: pressure, enthalpy, and specific humidities.
• pρ: pressure, density, and specific humidities.
• pθ_li: pressure, liquid-ice potential temperature, and specific humidities.
• ρθ_li: density, liquid-ice potential temperature, and specific humidities.

ρe <: IndepVars

Density, internal energy, and specific humidity

pe <: IndepVars

Pressure, internal energy, and specific humidity

ph <: IndepVars

Pressure, enthalpy, and specific humidity

pρ <: IndepVars

Pressure, density, and specific humidity

pθ_li <: IndepVars

Pressure, liquid-ice potential temperature, and specific humidity

ρθ_li <: IndepVars

Density, liquid-ice potential temperature, and specific humidity

DryAdiabaticProcess

A type for dispatching to isentropic (dry adiabatic) formulas.

gas_constant_air(param_set, q_tot=0, q_liq=0, q_ice=0)

The specific gas constant of moist air R_m.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• R_m: specific gas constant of moist air [J/(kg·K)]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression.

cp_m(param_set, q_tot=0, q_liq=0, q_ice=0)

The isobaric specific heat capacity of moist air cp_m.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• q_tot: total specific humidity of water [kg/kg]
• q_liq: specific humidity of liquid [kg/kg]
• q_ice: specific humidity of ice [kg/kg]

Returns

• cp_m: isobaric specific heat capacity [J/(kg·K)]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The specific heat capacities are assumed to be constant (calorically perfect air).

cv_m(param_set, q_tot=0, q_liq=0, q_ice=0)

The isochoric specific heat capacity of moist air cv_m.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• cv_m: isochoric specific heat capacity [J/(kg·K)]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The specific heat capacities are assumed to be constant (calorically perfect air).

soundspeed_air(param_set, T, q_tot=0, q_liq=0, q_ice=0)

The speed of sound in unstratified air.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• c: speed of sound [m/s]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The formula is c = √(γ R_m T) where γ = cp_m/cv_m.

vapor_specific_humidity(q_tot=0, q_liq=0, q_ice=0)

The vapor specific humidity.

Arguments

• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• q_vap: vapor specific humidity [kg/kg]

If the specific humidities are not given, the result is zero. The formula is q_vap = q_tot - q_liq - q_ice.

Note: callers must ensure q_liq + q_ice ≤ q_tot; no clamping is applied.

condensate_specific_humidity(q_liq=0, q_ice=0)

The condensate specific humidity.

Arguments

• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• q_cond: condensate specific humidity [kg/kg]

If the specific humidities are not given, the result is zero. The formula is q_cond = q_liq + q_ice.

vol_vapor_mixing_ratio(param_set, q_tot=0, q_liq=0, q_ice=0)

The molar (volume) mixing ratio of water vapor.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• r_v: molar mixing ratio of water vapor [mol/mol]

If the specific humidities are not given, the result is zero.

partial_pressure_dry(param_set, p, q_tot=0, q_liq=0, q_ice=0)

The partial pressure of dry air.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• p: air pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• p_d: partial pressure of dry air [Pa]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this equals the total pressure.

partial_pressure_vapor(param_set, p, q_tot=0, q_liq=0, q_ice=0)

The partial pressure of water vapor.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• p: air pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• p_v: partial pressure of water vapor [Pa]

If the specific humidities are not given, the partial pressure is zero.

vapor_pressure_deficit(param_set, T, p, q_tot, q_liq, q_ice, phase::Phase)

The vapor pressure deficit (saturation vapor pressure minus actual vapor pressure, truncated to be non-negative) over a specific phase (Liquid or Ice).

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: air temperature [K]
• p: air pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]
• phase: phase to compute saturation over (Liquid() or Ice())

Returns

• VPD: vapor pressure deficit [Pa]

vapor_pressure_deficit(param_set, T, p, q_tot=0, q_liq=0, q_ice=0)

The vapor pressure deficit (saturation vapor pressure minus actual vapor pressure, truncated to be non-negative) over liquid water for temperatures above freezing and over ice for temperatures below freezing.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: air temperature [K]
• p: air pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• VPD: vapor pressure deficit [Pa]

If the specific humidities are not given, the result is the saturation vapor pressure.

relative_humidity(param_set, T, p, q_tot=0, q_liq=0, q_ice=0)

The relative humidity (clipped between 0 and 1).

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• p: pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• RH: relative humidity [dimensionless], 0 ≤ RH ≤ 1

If q_liq and q_ice are zero (or not given), the relative humidity is computed relative to saturation over ice below freezing and over liquid above freezing. If condensate is present, the relative humidity is computed relative to saturation over a mixture of liquid and ice, with the liquid fraction given by the ratio q_liq / (q_liq + q_ice).

Note: relative_humidity uses saturation_vapor_pressure(param_set, T, q_liq, q_ice). In particular, for q_liq == q_ice == 0 it includes a small smooth transition around freezing (via liquid_fraction).

condensate_partition(param_set, T, ρ, q_tot)

Compute the equilibrium liquid and ice specific humidities from the condensate in thermodynamic equilibrium, returning a tuple (q_liq, q_ice).

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]

Returns

• (q_liq, q_ice): tuple of liquid and ice specific humidities [kg/kg]

The condensate is partitioned into liquid and ice using the temperature-dependent liquid fraction from liquid_fraction_ramp (a power-law interpolation between T_icenuc and T_freeze) and the saturation excess (see saturation_excess).

has_condensate(param_set, q_c)

Bool indicating if condensate exists, i.e., q_c > q_min.

We use the small threshold q_min rather than 0 to avoid division by zero in functions such as liquid_fraction and to robustly handle numerical noise.

liquid_fraction(param_set, T, q_liq, q_ice)

The fraction of condensate that is liquid.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• λ: liquid fraction [dimensionless], 0 ≤ λ ≤ 1

The liquid fraction is computed from q_liq and q_ice. If q_liq + q_ice exceeds a small threshold (see has_condensate), q_liq / (q_liq + q_ice) is returned. If there is effectively no condensate, a smooth temperature-dependent partitioning is used: a linear ramp from 0 to 1 over the 0.2 K interval [T_freeze - 0.2, T_freeze], so that the liquid fraction is exactly 1 at T_freeze.

liquid_fraction_ramp(param_set, T)

The liquid fraction computed from temperature using a power law interpolation between T_icenuc and T_freeze.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]

Returns

• λ: liquid fraction [dimensionless], 0 ≤ λ ≤ 1

Edge cases:

• For T > T_freeze, this returns 1; for T ≤ T_icenuc, it returns 0.

Reference

Kaul et al. (2015), "Sensitivities in large-eddy simulations of mixed-phase Arctic stratocumulus clouds using a simple microphysics approach," Monthly Weather Review, 143, 4393-4421, doi:10.1175/MWR-D-14-00319.1.

specific_humidity_to_mixing_ratio(q, q_tot)

Converts specific humidity to mixing ratio (total basis).

Arguments

• q: specific humidity of interest [kg/kg]
• q_tot: total specific humidity [kg/kg]

Returns

• r: mixing ratio [kg/kg]

Note that this function is singular when q_tot = 1. The formula is r = q / (1 - q_tot).

internal_energy(param_set, T, q_tot=0, q_liq=0, q_ice=0)

The internal energy per unit mass.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• e_int: specific internal energy [J/kg]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The internal energy is computed as a mass-weighted sum of the internal energies of each component (dry air, vapor, liquid, ice), referenced to T_0.

internal_energy_dry(param_set, T)

The dry air internal energy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]

Returns

• e_d: specific internal energy of dry air [J/kg]

internal_energy_vapor(param_set, T)

The water vapor internal energy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]

Returns

• e_v: specific internal energy of water vapor [J/kg]

internal_energy_liquid(param_set, T)

The liquid water internal energy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]

Returns

• e_l: specific internal energy of liquid water [J/kg]

internal_energy_ice(param_set, T)

The ice internal energy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]

Returns

• e_i: specific internal energy of ice [J/kg]

enthalpy(e_int, R_m, T)

The specific enthalpy.

Arguments

• e_int: internal specific energy [J/kg]
• R_m: gas constant of moist air [J/(kg·K)], see gas_constant_air
• T: air temperature [K]

Returns

• h: specific enthalpy [J/kg]

The enthalpy is computed as h = e_int + R_m T, which follows from h = e_int + p v with the ideal gas law p v = R_m T (specific volume v = 1/ρ).

enthalpy(param_set, T, q_tot=0, q_liq=0, q_ice=0)

The specific enthalpy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• h: specific enthalpy [J/kg]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The enthalpy is computed as a mass-weighted sum of the enthalpies of each component (dry air, vapor, liquid, ice).

enthalpy_dry(param_set, T)

The specific enthalpy of dry air.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]

Returns

• h_d: specific enthalpy of dry air [J/kg]

enthalpy_vapor(param_set, T)

The specific enthalpy of water vapor.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]

Returns

• h_v: specific enthalpy of water vapor [J/kg]

enthalpy_liquid(param_set, T)

The specific enthalpy of liquid.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]

Returns

• h_l: specific enthalpy of liquid [J/kg]

The specific enthalpy of liquid is equal to the internal energy of liquid because the specific volume of condensed water is neglected (i.e., p v_l ≈ 0).

enthalpy_ice(param_set, T)

The specific enthalpy of ice.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]

Returns

• h_i: specific enthalpy of ice [J/kg]

The specific enthalpy of ice is equal to the internal energy of ice because the specific volume of condensed water is neglected (i.e., p v_i ≈ 0).

total_energy(param_set, e_kin, e_pot, T, q_tot=0, q_liq=0, q_ice=0)

The total energy per unit mass.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• e_kin: kinetic energy per unit mass [J/kg]
• e_pot: gravitational potential energy (geopotential) per unit mass [J/kg]
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• e_tot: total specific energy [J/kg]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The total energy is e_tot = e_int + e_kin + e_pot.

total_enthalpy(param_set, e_tot, T, q_tot=0, q_liq=0, q_ice=0)

The total enthalpy per unit mass, h_tot = e_tot + R_m T.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• e_tot: total specific energy [J/kg], see total_energy
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• h_tot: total specific enthalpy [J/kg]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression.

dry_static_energy(param_set, T, e_pot)

The dry static energy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• e_pot: gravitational potential energy per unit mass (geopotential) [J/kg]

Returns

• s_d: dry static energy [J/kg]

The dry static energy is the sum of the dry enthalpy and geopotential: s_d = h_d + e_pot.

vapor_static_energy(param_set, T, e_pot)

The static energy (sensible heat only) of water vapor.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• e_pot: gravitational potential energy per unit mass (geopotential) [J/kg]

Returns

• s_v: vapor static energy [J/kg]

The formula is s_v = cp_v * (T - T_0) + e_pot, where T_0 is the reference temperature.

moist_static_energy(param_set, T, e_pot, q_tot=0, q_liq=0, q_ice=0)

The moist static energy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• e_pot: gravitational potential energy per unit mass (geopotential) [J/kg]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• s_m: moist static energy [J/kg]

If the specific humidities are not given, the result is for dry air. The moist static energy is the sum of the moist enthalpy and geopotential: s_m = h + e_pot, where h is computed from enthalpy.

virtual_dry_static_energy(param_set, T, e_pot, q_tot=0, q_liq=0, q_ice=0)

The virtual dry static energy.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• e_pot: gravitational potential energy per unit mass (geopotential) [J/kg]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• s_vd: virtual dry static energy [J/kg]

If the specific humidities are not given, the result is for dry air. The virtual dry static energy is s_vd = cp_d * (T_virt - T_0) + e_pot, where T_virt is the virtual temperature (see virtual_temperature) and T_0 is the reference temperature.

latent_heat_fusion(param_set, T)

The specific latent heat of fusion L_f.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]

Returns

• L_f: latent heat of fusion [J/kg]

Because the specific heats are assumed constant, the latent heats are linear functions of temperature by Kirchhoff's law.

latent_heat_sublim(param_set, T)

The specific latent heat of sublimation L_s.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]

Returns

• L_s: latent heat of sublimation [J/kg]

Because the specific heats are assumed constant, the latent heats are linear functions of temperature by Kirchhoff's law.

latent_heat_vapor(param_set, T)

The specific latent heat of vaporization L_v.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]

Returns

• L_v: latent heat of vaporization [J/kg]

Because the specific heats are assumed constant, the latent heats are linear functions of temperature by Kirchhoff's law.

latent_heat_mixed(param_set, T, λ)

The specific latent heat of a mixed phase, weighted by the liquid fraction λ.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: air temperature [K]
• λ: liquid fraction [dimensionless], 0 ≤ λ ≤ 1

Returns

• L: latent heat of the mixed phase [J/kg]

Computed as L = λ * L_v(T) + (1 - λ) * L_s(T).

latent_heat(param_set, T)
latent_heat(param_set, T, q_liq, q_ice)

The specific latent heat of phase change for a mixed-phase condensate.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• q_liq: (optional) liquid specific humidity [kg/kg]
• q_ice: (optional) ice specific humidity [kg/kg]

Returns

• L: latent heat [J/kg]

If q_liq and q_ice are provided, the liquid fraction is computed from these (see liquid_fraction). If there is no condensate, a temperature-dependent partitioning similar to liquid_fraction_ramp is used.

If q_liq and q_ice are not provided, the liquid fraction is computed from a temperature-dependent parameterization liquid_fraction_ramp.

See also latent_heat_mixed for the underlying computation.

humidity_weighted_latent_heat(param_set, q_liq=0, q_ice=0)

Specific-humidity weighted sum of latent heats of liquid and ice evaluated at reference temperature T_0.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• L_q: humidity-weighted latent heat [J/kg]

Computed as L_q = LH_v0 * q_liq + LH_s0 * q_ice. If q_liq and q_ice are not provided, this returns zero. Used in liquid-ice potential temperature computations.

air_temperature(param_set, e_int, q_tot=0, q_liq=0, q_ice=0)

The air temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• e_int: specific internal energy of moist air [J/kg]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• T: air temperature [K]

If the specific humidities are not given, the result is for dry air. This method inverts internal_energy by solving for T given e_int.

air_temperature(param_set, ::ph, h, q_tot=0, q_liq=0, q_ice=0)

The air temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• h: specific enthalpy of moist air [J/kg]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• T: air temperature [K]

If the specific humidities are not given, the result is for dry air. This method inverts enthalpy by solving for T given h.

air_temperature(param_set, ::pρ, p, ρ, q_tot=0, q_liq=0, q_ice=0)

The air temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• p: air pressure [Pa]
• ρ: air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• T: air temperature [K]

If the specific humidities are not given, the result is for dry air. This directly applies the ideal gas law: T = p / (R_m ρ).

air_temperature(
    param_set,
    ::pθ_li,
    p,
    θ_li,
    q_tot=0,
    q_liq=0,
    q_ice=0
)

The air temperature obtained by inverting the liquid-ice potential temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• p: pressure [Pa]
• θ_li: liquid-ice potential temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• T: air temperature [K]

If the specific humidities are not given, θ_li is assumed to be the dry-air potential temperature. This inverts liquid_ice_pottemp_given_pressure by solving for T.

air_temperature(param_set, ::ρθ_li, ρ, θ_li, q_tot=0, q_liq=0, q_ice=0)

The air temperature obtained by inverting the liquid-ice potential temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• ρ: (moist-)air density [kg/m³]
• θ_li: liquid-ice potential temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• T: air temperature [K]

If the specific humidities are not given, the results are for dry air.

Method

This function uses a second-order Taylor expansion to avoid an iterative inversion. The unsaturated temperature T_unsat is first computed assuming the ideal gas law with potential temperature, then latent heat corrections are applied.

air_temperature(param_set, ::DryAdiabaticProcess, p, θ)

The air temperature for an isentropic process.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• p: pressure [Pa]
• θ: potential temperature [K]

Returns

• T: air temperature [K]

The temperature is computed using the definition of the dry potential temperature: T = θ * (p/p₀)^κ_d, where p₀ is the reference pressure and κ_d = R_d/cp_d.

potential_temperature(param_set, T, ρ, q_tot=0, q_liq=0, q_ice=0)

The potential temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• θ: potential temperature [K]

If the specific humidities are not given, the result is for dry air.

Note: if any of q_tot, q_liq, or q_ice are nonzero, the Exner function uses the moist exponent R_m/cp_m (reducing to the dry exponent R_d/cp_d in the dry limit).

potential_temperature_given_pressure(param_set, T, p, q_tot=0, q_liq=0, q_ice=0)

The potential temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• p: pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• θ: potential temperature [K]

If the specific humidities are not given, the result is for dry air.

Note: if any of q_tot, q_liq, or q_ice are nonzero, the Exner function uses the moist exponent R_m/cp_m (reducing to the dry exponent R_d/cp_d in the dry limit).

liquid_ice_pottemp(param_set, T, ρ, q_tot=0, q_liq=0, q_ice=0)

The liquid-ice potential temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• θ_li: liquid-ice potential temperature [K]

If the specific humidities are not given, the result is for dry air.

Reference

See [9] and the Mathematical Formulation §12.4.

liquid_ice_pottemp_given_pressure(param_set, T, p, q_tot=0, q_liq=0, q_ice=0)

The liquid-ice potential temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• p: pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• θ_li: liquid-ice potential temperature [K]

If the specific humidities are not given, the result is for dry air. The latent heats of phase transitions are approximated as constants.

Reference

See [9] and the Mathematical Formulation §12.4.

virtual_pottemp(param_set, T, ρ, q_tot=0, q_liq=0, q_ice=0)

The virtual potential temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• θ_v: virtual potential temperature [K]

If the specific humidities are not given, the result is for dry air. The virtual potential temperature is defined as θ_v = θ * R_m / R_d.

virtual_temperature(param_set, T, q_tot=0, q_liq=0, q_ice=0)

The virtual temperature.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• T_v: virtual temperature [K]

If the specific humidities are not given, the result is for dry air. The virtual temperature is defined such that dry air at T_v has the same density as moist air at T, i.e., T_v = T (R_m / R_d).

air_pressure(param_set, T, ρ, q_tot=0, q_liq=0, q_ice=0)

The air pressure from the equation of state (ideal gas law).

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: air temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• p: air pressure [Pa]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The formula is p = ρ R_m T where R_m is the gas constant of moist air (see gas_constant_air).

air_pressure(param_set, ::DryAdiabaticProcess, T, T∞, p∞)

The air pressure for an isentropic process.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: current temperature [K]
• T∞: reference temperature [K]
• p∞: reference pressure [Pa]

Returns

• p: air pressure [Pa]

The pressure is computed using the isentropic relation: p = p∞ * (T/T∞)^(1/κ_d), where κ_d = R_d/cp_d is the ratio of the gas constant to the isobaric specific heat constant of dry air.

air_density(param_set, T, p, q_tot=0, q_liq=0, q_ice=0)

The (moist-)air density from the equation of state (ideal gas law).

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: air temperature [K]
• p: pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• ρ: (moist-)air density [kg/m³]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The formula is ρ = p / (R_m T) where R_m is the gas constant of moist air (see gas_constant_air).

exner(param_set, T, ρ, q_tot=0, q_liq=0, q_ice=0)

The Exner function Π = (p/p₀)^(R_m/cp_m).

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• Π: Exner function [dimensionless]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The pressure is computed internally from the equation of state (see air_pressure).

exner_given_pressure(param_set, p, q_tot=0, q_liq=0, q_ice=0)

The Exner function Π = (p/p₀)^(R_m/cp_m), where p₀ is the reference pressure.

Arguments

• param_set: thermodynamics parameter set, see the Thermodynamics for more details
• p: pressure [Pa]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• Π: Exner function [dimensionless]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The reference pressure p₀ is p_ref_theta from the parameter set.

air_pressure_given_θ(param_set, ::DryAdiabaticProcess, θ, Φ)

The air pressure for an isentropic process.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• θ: potential temperature [K]
• Φ: gravitational potential [J/kg]

Returns

• p: air pressure [Pa]

The pressure is computed using the hydrostatic balance and the definition of potential temperature for an isentropic process: p = p₀ * (1 - Φ/(θ * cp_d))^(1/κ_d), where p₀ is the reference pressure, cp_d is the isobaric specific heat capacity of dry air, and κ_d = R_d/cp_d.

saturation_vapor_pressure(param_set, T, ::Liquid)
saturation_vapor_pressure(param_set, T, ::Ice)

The saturation vapor pressure over a plane surface of condensate.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• Liquid() or Ice() to dispatch over the condensate type

Returns

• p_v^*: saturation vapor pressure [Pa]

The saturation vapor pressure is computed by integration of the Clausius-Clapeyron relation, assuming constant specific heat capacities in the so-called Rankine-Kirchhoff approximation.

saturation_vapor_pressure(param_set, T)
saturation_vapor_pressure(param_set, T, q_liq, q_ice)

The saturation vapor pressure over liquid, ice, or a mixture of liquid and ice.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• q_liq: (optional) liquid specific humidity [kg/kg]
• q_ice: (optional) ice specific humidity [kg/kg]

Returns

• p_v^*: saturation vapor pressure [Pa]

Important: two different liquid fraction parameterizations are used

The two-argument form saturation_vapor_pressure(param_set, T) uses liquid_fraction_ramp, a power-law interpolation between T_icenuc and T_freeze.

The four-argument form saturation_vapor_pressure(param_set, T, q_liq, q_ice) uses liquid_fraction, which returns q_liq / (q_liq + q_ice) when condensate is present, or a linear ramp over [T_freeze - 0.2 K, T_freeze] when there is no condensate.

As a result, saturation_vapor_pressure(param_set, T, 0, 0) and saturation_vapor_pressure(param_set, T) give the same value at T_freeze (both λ=1) but can differ away from freezing. Use the two-argument form when no condensate information is available; use the four-argument form when condensate fractions are known.

q_vap_saturation(param_set, T, ρ)
q_vap_saturation(param_set, T, ρ, q_liq, q_ice)

The saturation specific humidity.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• ρ: air density [kg/m³]
• (optional) q_liq: liquid specific humidity [kg/kg]
• (optional) q_ice: ice specific humidity [kg/kg]

Returns

• q_v^*: saturation specific humidity [kg/kg]

If q_liq and q_ice are provided, the saturation specific humidity is that over a mixture of liquid and ice, computed in a thermodynamically consistent way from the weighted sum of the latent heats of the respective phase transitions. That is, the saturation vapor pressure and from it the saturation specific humidity are computed from a weighted mean of the latent heats of vaporization and sublimation, with the weights given by the liquid fraction (see liquid_fraction). If q_liq and q_ice are 0, the saturation specific humidity is that over liquid above freezing and over ice below freezing.

Otherwise, the fraction of liquid is given by the temperature dependent liquid_fraction_ramp(param_set, T).

Reference

Thermodynamically consistent mixed-phase saturation following [7].

q_vap_saturation(param_set, T, ρ, phase::Phase)

The saturation specific humidity.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• phase: the phase to compute saturation over (either Liquid() or Ice())

Returns

• q_v^*: saturation specific humidity [kg/kg]

∂q_vap_sat_∂T(param_set, T, ρ)
∂q_vap_sat_∂T(param_set, T, ρ, q_liq, q_ice)
∂q_vap_sat_∂T(param_set, T, ρ, phase::Phase)

Derivative of saturation vapor specific humidity with respect to temperature.

Arguments

• param_set: Parameter set.
• T: Temperature [K].
• ρ: Density [kg/m³].
• q_liq: (Optional) Liquid specific humidity [kg/kg].
• q_ice: (Optional) Ice specific humidity [kg/kg].
• phase: (Optional) Phase (Liquid() or Ice()).

Returns

• ∂q_vap_sat / ∂T: Derivative of saturation specific humidity with respect to temperature [kg/kg/K].

Computed as

\[\frac{\partial q_v^*}{\partial T} = q_v^* \left( \frac{L}{R_v T^2} - \frac{1}{T} \right),\]

which follows from the definition of saturation specific humidity and the ideal gas law,

\[q_v^* = \frac{p_v^*}{\rho R_v T},\]

(where $q_v^*$ is q_vap_sat and $p_v^*$ is saturation vapor pressure) and the Clausius-Clapeyron relation

\[\frac{\partial \ln p_v^*}{\partial T} = \frac{L}{R_v T^2}\]

by differentiation with respect to $T$ while holding $\rho$ constant.

If q_liq and q_ice are provided, the saturation specific humidity derivative is computed assuming a phase mixture defined by the liquid fraction (see liquid_fraction).

If phase is provided, the derivative is computed for saturation over that specific phase.

If neither are provided, the saturation specific humidity derivative is computed assuming phase equilibrium, where the liquid fraction is determined by a temperature-dependent parameterization (see liquid_fraction_ramp).

supersaturation(param_set, q_vap, ρ, T[, phase::Phase = Liquid()])

The supersaturation over water or ice.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• q_vap: vapor specific humidity [kg/kg]
• ρ: air density [kg/m³]
• T: air temperature [K]
• phase: (optional) liquid or ice phase to dispatch over (default: Liquid())

Returns

• S: supersaturation [dimensionless], defined as p_v/p_v^* - 1

supersaturation(param_set, q_vap, ρ, T, p_v_sat)

The supersaturation given the saturation vapor pressure.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• q_vap: vapor specific humidity [kg/kg]
• ρ: air density [kg/m³]
• T: temperature [K]
• p_v_sat: saturation vapor pressure [Pa]

Returns

• S: supersaturation [dimensionless], defined as p_v/p_v_sat - 1

saturation_excess(param_set, T, ρ, q_tot)
saturation_excess(param_set, T, ρ, q_tot, q_liq, q_ice)

The saturation excess in equilibrium.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• (optional) q_liq: liquid specific humidity [kg/kg]
• (optional) q_ice: ice specific humidity [kg/kg]

Returns

• q_ex: saturation excess [kg/kg]

The saturation excess is the difference between the total specific humidity q_tot and the saturation specific humidity in equilibrium, and it is defined to be nonzero only if this difference is positive: q_ex = max(0, q_tot - q_v^*).

saturation_excess(param_set, T, ρ, q_tot, p_vap_sat)

The saturation excess given the saturation vapor pressure.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• ρ: air density [kg/m³]
• q_tot: total specific humidity [kg/kg]
• p_vap_sat: saturation vapor pressure [Pa]

Returns

• q_ex: saturation excess [kg/kg]

The saturation excess is the difference between the total specific humidity q_tot and the saturation specific humidity, and it is defined to be nonzero only if this difference is positive: q_ex = max(0, q_tot - q_v^*).

q_vap_from_RH(param_set, p, T, RH, phase)

Compute the vapor specific humidity from the relative humidity.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• p: pressure [Pa]
• T: temperature [K]
• RH: relative humidity [dimensionless], 0 ≤ RH ≤ 1
• phase: the phase to compute saturation over (either Liquid() or Ice())

Returns

• q_vap: vapor specific humidity [kg/kg]

Note: this function assumes all water is in vapor form (no liquid or ice condensate).

q_vap_from_p_vap(param_set, T, ρ, p_v)

Compute the vapor specific humidity from the vapor partial pressure.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: air temperature [K]
• ρ: (moist-)air density [kg/m³]
• p_v: partial pressure of water vapor [Pa]

Returns

• q_vap: vapor specific humidity [kg/kg]

q_vap_saturation_from_pressure(param_set, q_tot, p, T)
q_vap_saturation_from_pressure(param_set, q_tot, p, T, q_liq, q_ice)

The saturation specific humidity.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• q_tot: total water specific humidity [kg/kg]
• p: air pressure [Pa]
• T: air temperature [K]
• (optional) q_liq: liquid specific humidity [kg/kg]
• (optional) q_ice: ice specific humidity [kg/kg]

Returns

• q_v^*: saturation specific humidity [kg/kg]

If q_liq and q_ice are provided, the saturation vapor pressure is computed from a weighted mean of the latent heats of vaporization and sublimation, with the weights given by the liquid fraction (see liquid_fraction). If q_liq and q_ice are 0, the saturation vapor pressure is that over liquid above freezing and over ice below freezing.

Otherwise, the liquid fraction is computed from a temperature dependent parameterization liquid_fraction_ramp(param_set, T).

The saturation specific humidity is computed as: $q_v^* = (R_d / R_v) * (1 - q_{tot}) * p_v^* / (p - p_v^*)$ where p_v^* is the saturation vapor pressure.

Edge case: this expression assumes p > p_v^*(T); if p ≤ p_v^*(T) the denominator changes sign.

saturation_adjustment(
    ::Type{M},  # RS.RootSolvingMethod type
    param_set,
    ::ρe,    ρ, e_int, q_tot, maxiter, tol, [T_guess], [forced_fixed_iters]
)
saturation_adjustment(
    ::Type{M}, param_set,
    ::pe,    p, e_int, q_tot, maxiter, tol, [T_guess], [forced_fixed_iters]
)
saturation_adjustment(
    ::Type{M}, param_set,
    ::ph,    p, h,     q_tot, maxiter, tol, [T_guess], [forced_fixed_iters]
)
saturation_adjustment(
    ::Type{M}, param_set,
    ::pρ,    p, ρ,     q_tot, maxiter, tol, [T_guess], [forced_fixed_iters]
)
saturation_adjustment(
    ::Type{M}, param_set,
    ::pθ_li, p, θ_li,  q_tot, maxiter, tol, [T_guess], [forced_fixed_iters]
)
saturation_adjustment(
    ::Type{M}, param_set,
    ::ρθ_li, ρ, θ_li,  q_tot, maxiter, tol, [T_guess], [forced_fixed_iters]
)

Solve for the temperature at which the given state is in phase equilibrium, and return that temperature together with the resulting condensate.

Arguments

• M: root-solving method from RootSolvers, e.g. RS.NewtonsMethod or RS.SecantMethod. RS.NewtonsMethod is recommended: analytic derivatives are available for every formulation.
• param_set: thermodynamics parameter set, see Thermodynamics
• indep_vars: an IndepVars singleton naming the two independent variables
• var₁, var₂: those two variables, in the order the type name spells them. Units follow the variable: ρ [kg/m³], p [Pa], e_int and h [J/kg], θ_li [K].
• q_tot: total specific humidity [kg/kg]
• maxiter: maximum number of solver iterations
• tol: relative tolerance on the temperature, or an RS.RelativeSolutionTolerance
• T_guess: optional initial temperature guess [K]. Defaults to nothing.
• forced_fixed_iters: run maxiter iterations without testing convergence, avoiding branch divergence on GPUs. T_guess and tol are ignored when true. Defaults to false.

Returns

• NamedTuple (; T, q_liq, q_ice, converged):
  • T: temperature [K]
  • q_liq: liquid specific humidity [kg/kg]
  • q_ice: ice specific humidity [kg/kg]
  • converged: whether the solver converged

Notes

• The state is first tested against the temperature it would have with no condensate; if that temperature is already subsaturated it is returned directly, since it solves the problem exactly. Otherwise the equilibrium condition is solved by root-finding and the condensate follows from condensate_partition.
• GPU broadcasting: pass forced_fixed_iters as a positional Bool.

Examples

import RootSolvers as RS
import Thermodynamics as TD
using ClimaParams

param_set = TD.Parameters.ThermodynamicsParameters(Float64)

# A cloudy state near 290 K and 850 hPa. Internal energy is measured relative to the
# reference temperature T_0, so it is negative here.
ρ, e_int, q_tot = 1.0175, -30923.0, 0.019
sol = TD.saturation_adjustment(
    RS.NewtonsMethod, param_set, TD.ρe(), ρ, e_int, q_tot, 20, 1e-4,
)
sol.T, sol.q_liq, sol.q_ice, sol.converged   # ≈ (290.0, 0.0049, 0.0, true)

See also the convenience methods below, which pick GPU-friendly defaults for you.

saturation_adjustment(param_set, ::ρe,    ρ, e_int, q_tot; maxiter = 2)
saturation_adjustment(param_set, ::pe,    p, e_int, q_tot; maxiter = 2)
saturation_adjustment(param_set, ::ph,    p, h,     q_tot; maxiter = 2)
saturation_adjustment(param_set, ::pρ,    p, ρ,     q_tot; maxiter = 2)
saturation_adjustment(param_set, ::pθ_li, p, θ_li,  q_tot; maxiter = 2)
saturation_adjustment(param_set, ::ρθ_li, ρ, θ_li,  q_tot; maxiter = 2)

Convenience methods with GPU-friendly defaults.

Runs a fixed number of safeguarded Newton iterations (see saturation_adjustment_fixed_iters) rather than a convergence-tested solve, so every GPU lane executes the same instructions.

The default maxiter = 2 is chosen for how saturation adjustment is used in a time-stepping model, where it is applied every step to a state that was near equilibrium at the previous one. What sets the iteration count is the gap between the unsaturated first guess and the solution, i.e. the warming from condensing the excess vapor. Across the tested profiles that gap never exceeds 4 K, and two iterations hold the temperature error below 2e-3 K. The error grows for states quenched from much further out of equilibrium — up to a few tenths of a K around a 10 K gap and of order 10 K beyond a 20 K gap, depending on the formulation — which is where a larger maxiter is worth passing.

Overshooting the iteration budget is no longer dangerous: the error is a truncation error that decreases monotonically with maxiter, not a stalled iteration. Earlier versions could settle into a limit cycle and return an answer tens of K off no matter how many iterations were allowed.

For more control over solver parameters, use the full signature with an explicit method type.

Returns

• NamedTuple (; T, q_liq, q_ice). Use the full signature if you need the converged flag.

Examples

import Thermodynamics as TD
using ClimaParams

param_set = TD.Parameters.ThermodynamicsParameters(Float64)
sol = TD.saturation_adjustment(param_set, TD.ρe(), 1.0175, -30923.0, 0.019)
sol.T, sol.q_liq, sol.q_ice   # ≈ (290.4, 0.0046, 0.0)
# This state condenses ~5 g/kg in one call (a ~13 K gap), so two iterations leave ~0.4 K;
# the converged answer is 290.0 K. See the accuracy discussion above.

internal_energy_sat(param_set, T, ρ, q_tot)

The internal energy per unit mass in thermodynamic equilibrium at saturation.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]

Returns

• e_int: specific internal energy [J/kg]

The phase partition into liquid and ice is computed internally from q_tot using the temperature-dependent liquid fraction (see liquid_fraction_ramp) and saturation excess (see saturation_excess).

enthalpy_sat(param_set, T, ρ, q_tot)

The specific enthalpy in thermodynamic equilibrium at saturation.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• T: temperature [K]
• ρ: (moist-)air density [kg/m³]
• q_tot: total specific humidity [kg/kg]

Returns

• h: specific enthalpy [J/kg]

The phase partition into liquid and ice is computed internally from q_tot using the temperature-dependent liquid fraction (see liquid_fraction_ramp) and saturation excess (see saturation_excess).

∂e_int_∂T_sat_ρ(param_set, T, ρ, q_tot, unsat_branch = Val(true))

Derivative of internal_energy_sat with respect to temperature at fixed density.

Uses ∂q_vap_sat/∂T|_ρ from Clausius-Clapeyron, which includes the -1/T term from the density dependence of saturation vapor pressure.

Passing Val(false) for unsat_branch returns the saturated form even below saturation; see _select_sat_branch.

∂e_int_∂T_sat_p(param_set, T, p, q_tot)

Derivative of internal_energy_sat with respect to temperature at fixed pressure.

Uses ∂q_vap_sat/∂T|_p = q_vap_sat * L / (R_v T²) (Clausius-Clapeyron at constant pressure), which differs from the constant-density form used in ∂e_int_∂T_sat_ρ by an additional q_vap_sat / T term.

∂θ_li_∂T_sat_ρ(param_set, T, ρ, q_tot)

Derivative of liquid_ice_pottemp (at saturation equilibrium) with respect to temperature at fixed density.

Structured like ∂θ_li_∂T_sat_p, but accounts for p = ρ R_m T varying with T. This introduces two corrections to ∂θ/∂T:

• an extra -α/T from ∂ln(p)/∂T|_ρ = 1/T + ∂R_m/∂T / R_m, and
• an extra -∂R_m/∂T / cp_m from the same.

The ∂q_vap_sat/∂T also differs: at fixed ρ it carries an extra -q_vap_sat / T relative to the fixed-p Clausius–Clapeyron form.

∂θ_li_∂T_sat_p(param_set, T, p, q_tot)

Derivative of liquid_ice_pottemp_given_pressure (at saturation equilibrium) with respect to temperature at fixed pressure.

Uses the product rule on θ_li = θ * (1 - L_c / (cp_m T)), differentiating each factor through the T-dependent phase partition.

∂p_∂T_sat_ρ(param_set, T, ρ, q_tot)

Derivative of air_pressure (at saturation equilibrium) with respect to temperature at fixed density.

From p = ρ R_m T and ∂R_m/∂T = R_v · ∂q_vap_sat/∂T:

∂p/∂T|_ρ = ρ · (R_m + T · ∂R_m/∂T)

Unsaturated: R_m is constant, so ∂p/∂T = ρ R_m.

∂h_∂T_sat_p(param_set, T, p, q_tot)

Derivative of enthalpy_sat with respect to temperature at fixed pressure.

Structured identically to ∂e_int_∂T_sat_p but with component enthalpies (cp_m instead of cv_m, enthalpy_vapor instead of internal_energy_vapor, etc.).

entropy(param_set, p, T, q_tot=0, q_liq=0, q_ice=0)

The specific entropy in thermodynamic equilibrium.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• p: pressure [Pa]
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• s: specific entropy [J/(kg·K)]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), this reduces to the dry-air expression. The entropy is computed as a mass-weighted sum of the entropies of each component (dry air, vapor, liquid, ice).

Reference

Mass-weighted moist-air entropy formulation following [7].

entropy_dry(param_set, p, T, q_tot=0, q_liq=0, q_ice=0)

The specific entropy of dry air at its partial pressure.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• p: total air pressure [Pa]
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• s_d: specific entropy of dry air [J/(kg·K)]

In the dry limit (q_tot = q_liq = q_ice = 0, the default), the dry-air partial pressure equals the total pressure. Note: entropy_dry diverges logarithmically as q_tot → 1 (since p_d → 0). See also the analogous warning in entropy_vapor.

The entropy reference pressure is MSLP, which is not the reference pressure p_ref_theta used by exner and the potential temperatures. The two reference states are independent and need not agree; entropies from this function are therefore not directly comparable to potential temperatures without accounting for the offset.

entropy_vapor(param_set, p, T, q_tot=0, q_liq=0, q_ice=0)

The specific entropy of water vapor at its partial pressure.

Arguments

• param_set: thermodynamics parameter set, see Thermodynamics
• p: total air pressure [Pa]
• T: temperature [K]
• q_tot: total specific humidity [kg/kg]
• q_liq: liquid specific humidity [kg/kg]
• q_ice: ice specific humidity [kg/kg]

Returns

• s_v: specific entropy of water vapor [J/(kg·K)]

Note: the entropy of water vapor diverges logarithmically as q_tot → 0 (since p_v → 0). The analogous divergence occurs in entropy_dry as q_tot → 1 (since p_d → 0). A small numerical regularization ϵ_numerics(FT) is added to the pressure before taking the logarithm to avoid floating-point exceptions, but results should not be trusted in these limits.

TemperatureProfile

Abstract type for temperature or virtual temperature reference profiles that can be used in atmosphere models.

Instances of this type are required to be callable objects with the following signature

T,p = (::TemperatureProfile)(param_set::APS, z::FT) where {FT}

where T is the temperature or virtual temperature (K), and p is the pressure (Pa).

DecayingTemperatureProfile{FT} <: TemperatureProfile{FT}

Virtual temperature profile that decays smoothly with height z from T_virt_surf to T_min_ref over height scale H_t (default: density scale height with T_virt_surf).

\[T_{\text{v}}(z) = \max(T_{\text{v, sfc}} − (T_{\text{v, sfc}} - T_{\text{v, min}}) \tanh(z/H_{\text{t}}), T_{\text{v, min}})\]

Fields

• T_virt_surf: Virtual temperature at surface (K)
• T_min_ref: Minimum virtual temperature at the top of the atmosphere (K)
• H_t: Height scale over which virtual temperature drops (m)

IsothermalProfile(param_set, T_virt)
IsothermalProfile(param_set, ::Type{FT<:Real})

Uniform virtual temperature profile implemented as a special case of DecayingTemperatureProfile.

DryAdiabaticProfile{FT} <: TemperatureProfile{FT}

Temperature profile with uniform potential temperature θ up to the height where a minimum temperature is reached.

Fields

• T_surface: Surface temperature (K)
• T_min_ref: Minimum temperature (K)

DataCollection

This module is designed to help judge the accuracy and performance for a particular formulation, tolerance, and or solver configuration, by providing tools to collect various statistics when saturation_adjustment is called.

Example:

import Thermodynamics as TD
import RootSolvers as RS

function do_work()
    # Calls TD.saturation_adjustment()..., possibly many times
end

TD.solution_type() = RS.VerboseSolution()
do_work()
TD.DataCollection.print_summary()