Column Schemes

struct RadiativeTransferColumn{NF, LW, SW, PC, TC, V<:AbstractArray{NF, 1}, G<:(ColumnGrid{NF, V} where V<:AbstractArray{NF, 1}), AP<:(AtmosphereProfile{NF, V} where V<:AbstractArray{NF, 1})}

Bundles everything a single-column radiation calculation needs: grid, atmosphere profile, surface state, longwave and shortwave schemes, the thermodynamic and physical constants, and pre-allocated buffers for the temperature tendency, layer transmissivity, and diagnostic fluxes.

The high-level API is solve_longwave! and solve_shortwave! on a RadiativeTransferColumn:

column = RadiativeTransferColumn(; grid, profile, surface)
solve_longwave!(column)
solve_shortwave!(column)
@show column.longwave_diagnostics.outgoing_longwave

Fields:

• grid: Column vertical grid (sigma coordinates)
• profile: Atmosphere profile (temperature, humidity, geopotential, surface pressure, rain rate)
• surface: Lower boundary state (SST/LST, albedos, emissivities, cos-zenith)
• longwave_scheme: Longwave scheme, e.g. AnalyticBandLongwave
• shortwave_scheme: Shortwave scheme, e.g. OneBandShortwave or TransparentShortwave
• physical_constants: Physical constants (gravity, heat capacity, Stefan–Boltzmann, solar constant)
• thermodynamic_constants: Thermodynamic constants (Clausius–Clapeyron parameters for saturation humidity)
• temperature_tendency: Per-layer temperature tendency written by solve_longwave! / solve_shortwave!
• transmissivity_scratch: Per-layer scratch for the shortwave transmissivity
• longwave_diagnostics: Scalar longwave diagnostics (OLR, surface up/down, ocean/land split)
• shortwave_diagnostics: Scalar shortwave diagnostics (TOA up, surface up/down, albedo, clouds)

reset!(
    column::RadiativeTransferColumn
) -> RadiativeTransferColumn

Zero the temperature tendency and scalar diagnostics on column so a fresh solve_longwave! / solve_shortwave! doesn't accumulate onto stale values.

radiative_heating!(column::RadiativeTransferColumn; reset=true, longwave=true, shortwave=true)

High-level analytic-band column update. This is a convenience wrapper around the component calls solve_longwave! and solve_shortwave!; host models can keep using those lower-level calls directly when they own their own vertical integrals or tendency insertion.

radiation_workspace(model, atmosphere; backend=nothing)

Construct reusable storage for repeated runtime calls. Host integrations may also pass their own arrays/views directly to component methods.

radiation_workspace(column::RadiativeTransferColumn)

The existing single-column object is already a reusable workspace: it owns the temperature-tendency vector, shortwave transmissivity scratch, and diagnostic objects used by the analytic-band paths.

struct AtmosphereProfile{NF, V<:AbstractArray{NF, 1}}

Column thermodynamic profile and lower boundary quantities read by the radiation solvers.

The arrays are indexed top-down: k = 1 is the top of the atmosphere, k = Nz is the bottom (surface-adjacent) layer.

struct ColumnGrid{NF, V<:AbstractArray{NF, 1}}

Vertical geometry for a single column expressed in sigma-pressure coordinates.

σ_full has length Nz and gives the midpoint of each layer. σ_half has length Nz + 1 and gives the layer interfaces. σ_thick = diff(σ_half) has length Nz.

struct SurfaceState{NF}

Lower-boundary properties needed by the column radiation solvers. NaN in sea_surface_temperature or land_surface_temperature means the column has no ocean / no land respectively (the corresponding contribution is skipped).

struct PhysicalConstants{NF}

Physical constants consumed by the column radiation solvers, the staged runtime (a ColumnAtmosphere carries one) and the RRTMGP adapter. Hosts pass their own values through the keyword constructor; the fields and their Earth defaults are

• gravity — g, 9.80665 m s⁻²
• heat_capacity — cᵖ, isobaric specific heat of dry air, 1004.64 J kg⁻¹ K⁻¹
• stefan_boltzmann — σ, 5.670374419e-8 W m⁻² K⁻⁴
• solar_constant — S₀, 1361 W m⁻²
• dry_air_molar_mass — mᵈ, 0.0289647 kg mol⁻¹ (the ecCKD/CKDMIP value)
• water_molar_mass — mᵛ, 0.01801528 kg mol⁻¹
• dry_air_gas_constant — Rᵈ, 287.05 J kg⁻¹ K⁻¹
• universal_gas_constant — ℛ, 8.31446261815324 J mol⁻¹ K⁻¹
• avogadro_number — Nᴬ, 6.02214076e23 mol⁻¹

The column schemes read constants duck-typed, so any object with the properties they use works in their place.

struct ThermodynamicConstants{NF}

Thermodynamic constants needed for saturation-humidity calculations used by the diagnostic cloud scheme.

default_earth_constants(
    _::Type{NF}
) -> NamedTuple{(:physical, :thermodynamic), <:Tuple{PhysicalConstants, ThermodynamicConstants}}

Sensible Earth defaults for the full set of physical constants needed by the shortwave solver (constants + thermodynamic constants).

mutable struct LongwaveDiagnostics{NF}

Column longwave diagnostic outputs written by solve_longwave!.

mutable struct ShortwaveDiagnostics{NF}

Column shortwave diagnostic outputs written by solve_shortwave!.

AbstractLongwaveScheme <: AbstractRadiationScheme

A column longwave radiative transfer scheme. Concrete subtypes implement solve_longwave!.

struct AnalyticBandLongwave{NF} <: AbstractLongwaveScheme

Williams (2026) "Simple Spectral Model" (SSM) for clear-sky longwave radiative transfer.

The scheme solves Schwarzschild's two-stream equations

dℐꜛ/dτ = ℐꜛ − πB(T)
dℐꜜ/dτ = πB(T) − ℐꜜ

for the upwelling (ℐꜛˡʷ) and downwelling (ℐꜜˡʷ) spectral longwave fluxes at each of Nwavenumbers evenly spaced wavenumbers between wavenumber_min and wavenumber_max, with analytic mass absorption coefficients for H₂O line (rotation + vibration–rotation + combination bands, water_vapor_line_absorption_reference), a two-band H₂O continuum (water_vapor_continuum_absorption_reference) and a Lorentzian CO₂ 15 μm bending mode (carbon_dioxide_absorption_reference). All reference constants are at (Tref, pref, RH_ref) = (260 K, 500 hPa, 100 %).

Fields and defaults follow Williams (2026), Table 1.

References:

• Williams (2026), J. Adv. Model. Earth Syst., doi:10.1029/2025MS005405.
• Armstrong (1968), doi:10.1016/0022-4073(68)90052-6 (diffusivity factor D).
• Mlawer et al. (1997), doi:10.1029/97JD00237 (continuum temperature scaling).

Fields:

• Nwavenumbers: Number of evenly spaced wavenumber quadrature points
• wavenumber_min: Minimum wavenumber of the spectral integration range [cm⁻¹]
• wavenumber_max: Maximum wavenumber of the spectral integration range [cm⁻¹]
• κ_rot: Peak absorption of the pure-rotation band κ_rot [m² kg⁻¹]
• l_rot: e-folding decay length of the rotation band l_rot [cm⁻¹]
• κ_vr: Peak absorption of the vibration–rotation band κ_vr [m² kg⁻¹]
• l_vr1: e-folding length of vibration–rotation band (low-ν side) l_vr1 [cm⁻¹]
• l_vr2: e-folding length of vibration–rotation band (high-ν side) l_vr2 [cm⁻¹]
• κ_cnt1: Continuum absorption below 1700 cm⁻¹ κ_cnt1 [m² kg⁻¹]
• κ_cnt2: Continuum absorption above 1700 cm⁻¹ κ_cnt2 [m² kg⁻¹]
• κ_CO₂: Peak absorption of the CO₂ 15 μm band κ_CO₂ [m² kg⁻¹]
• l_CO₂: e-folding half-width of CO₂ band l_CO₂ [cm⁻¹]
• ν̃_CO₂: Centre wavenumber of CO₂ bending mode ν̃_CO₂ [cm⁻¹]
• diffusivity: Two-stream diffusivity factor D (Armstrong 1968)
• p_ref: Reference pressure for pressure broadening [Pa]
• T_ref: Reference temperature for absorption coefficient fits [K]
• pv_ref: Reference saturation water-vapor pressure at T_ref [Pa]
• σ_cont: Temperature-scaling exponent for the continuum (Mlawer et al. 1997) [K⁻¹]
• water_vapor_molar_mass_ratio: Water-to-dry-air molar mass ratio mᵛ/mᵈ of the vapor partial pressure
• carbon_dioxide_molar_mass_ratio: CO₂-to-dry-air molar mass ratio converting ppmv to a mass mixing ratio

solve_longwave!(
    temperature_tendency::AbstractVector,
    diagnostics::LongwaveDiagnostics{NF},
    scheme::AnalyticBandLongwave{NF},
    profile::AtmosphereProfile{NF, V} where V<:AbstractArray{NF, 1},
    geometry::ColumnGrid,
    surface::SurfaceState,
    constants
)

Column longwave radiative transfer for the Williams (2026) Simple Spectral Model. Tendencies are accumulated into temperature_tendency with +=/-=; diagnostic fluxes are written into diagnostics.

Sign convention: temperature tendency has units [K s⁻¹]; positive OLR, positive downward surface flux, positive upward surface flux.

solve_longwave!(
    column::RadiativeTransferColumn
) -> RadiativeTransferColumn

Column longwave radiative transfer using the scheme stored on column. Accumulates into column.temperature_tendency (call reset!(column) first if you want a clean slate) and writes scalars into column.longwave_diagnostics.

planck_wavenumber(T, ν̃)

Spectral Planck radiance at temperature T [K] and wavenumber ν̃ [cm⁻¹], in units of W m⁻² sr⁻¹ (cm⁻¹)⁻¹.

The two-stream hemispherical flux source term is π × planck_wavenumber(T, ν̃). Integrating π × planck_wavenumber(T, ν̃) over all wavenumbers recovers the Stefan–Boltzmann law σ T⁴.

The constants are the module's PLANCK_CONSTANT, SPEED_OF_LIGHT and BOLTZMANN_CONSTANT (CODATA 2018).

water_vapor_line_absorption_reference(ν̃, scheme::AnalyticBandLongwave) -> κ  [m² kg⁻¹]

Reference H₂O line mass absorption coefficient. Piecewise-exponential fit to the pure-rotation (200–1000 cm⁻¹), vibration–rotation (1000–1700 cm⁻¹) and combination bands (1700–2500 cm⁻¹).

Reference: Williams (2026), Eq. 4 and Table 1.

water_vapor_continuum_absorption_reference(ν̃, scheme::AnalyticBandLongwave) -> κ  [m² kg⁻¹]

Reference H₂O continuum absorption. Two gray values split at 1700 cm⁻¹ (stronger in the main atmospheric window below).

Reference: Williams (2026), Eq. 6.

carbon_dioxide_absorption_reference(ν̃, scheme::AnalyticBandLongwave) -> κ  [m² kg⁻¹]

Reference CO₂ absorption coefficient. A two-sided exponential (Laplace-shaped) wing centred on the 15 μm bending mode at ν̃_CO₂ ≈ 667 cm⁻¹, active only in [500, 850] cm⁻¹.

Reference: Williams (2026), Eq. 5.

williams_optical_depth_increment(
    k::Integer,
    ν̃,
    CO₂,
    temperature::AbstractVector,
    humidity::AbstractVector,
    surface_pressure::Real,
    geometry::ColumnGrid,
    scheme,
    gravity::Real
) -> Any

Optical depth increment through layer k at wavenumber ν̃. Combines H₂O line absorption (Williams 2026, Eq. 7), H₂O continuum (Eq. 8) and CO₂ (Eq. 9). The result already includes the two-stream diffusivity factor D ≈ 1.5 (Armstrong 1968).

temperature and humidity are length-Nz column vectors; geometry is a ColumnGrid. The molar mass ratios of the vapor partial pressure and of the CO₂ mass mixing ratio are the scheme's water_vapor_molar_mass_ratio and carbon_dioxide_molar_mass_ratio.

AbstractShortwaveScheme <: AbstractRadiationScheme

A column shortwave radiative transfer scheme. Concrete subtypes implement solve_shortwave!.

struct TransparentShortwave <: AbstractShortwaveScheme

Zero-atmosphere shortwave: TOA insolation reaches the surface unattenuated and is reflected by the surface albedo. Temperature tendencies are zero. Useful as a baseline and for tests of the surface energy budget.

struct OneBandShortwave{C<:AbstractShortwaveClouds, T<:AbstractShortwaveTransmissivity, R<:OneBandShortwaveRadiativeTransfer} <: AbstractShortwaveScheme

Composite one-band shortwave scheme: a cloud diagnosis, a layer-transmissivity model and a radiative-transfer solver are combined into a single scheme that can be passed to solve_shortwave!.

Fields:

• clouds: Cloud diagnosis, an AbstractShortwaveClouds
• transmissivity: Layer-transmissivity model, an AbstractShortwaveTransmissivity
• radiative_transfer: Radiative-transfer solver, an OneBandShortwaveRadiativeTransfer

OneBandGreyShortwave(
    ::Type{NF};
    clouds,
    transmissivity,
    radiative_transfer
) -> OneBandShortwave{NoClouds, T, R} where {T<:ConstantShortwaveTransmissivity, R<:(OneBandShortwaveRadiativeTransfer{_A, F} where {_A, F<:(NumericalRadiation.var"#34#35")})}

Convenience constructor with a dry-atmosphere default (no clouds, constant column transmissivity).

struct OneBandShortwaveRadiativeTransfer{NF, F} <: AbstractShortwaveScheme

SPEEDY-style one-band shortwave radiative transfer (Kucharski, Molteni & Bracco, 2006, Appendix). Computes cloud albedo reflection at the diagnosed cloud top, ozone absorption in the stratosphere, layer-by-layer transmission with the configured AbstractShortwaveTransmissivity, stratocumulus reflection just above the surface, and surface-albedo reflection.

Fields:

• ozone_absorption: Total ozone absorption as a fraction of incoming TOA flux (default NF(0.01))
• ozone_distribution: Ozone vertical distribution ζ(sigma_level) → weight over the sigma coordinate, normalised so ∫ ζ dσ = 1 (default default_ozone_distribution(NF))

AbstractShortwaveTransmissivity <: AbstractRadiationScheme

A layer-transmissivity model used by shortwave radiative-transfer solvers.

struct ConstantShortwaveTransmissivity{NF} <: AbstractShortwaveTransmissivity

Constant atmospheric transmissivity, distributed across layers proportional to their pressure thickness.

Fields:

• transmissivity: Column-integrated atmospheric transmissivity (0 .. 1) (default NF(0.85))

struct BackgroundShortwaveTransmissivity{NF} <: AbstractShortwaveTransmissivity

SPEEDY-style background shortwave transmissivity (Kucharski, Molteni & Bracco, 2006; cf. Fortran SPEEDY absdry, absaer, abswv1/abswv2, abscl1/abscl2, azen, nzen). The layer optical depth sums contributions from dry air, aerosols (∝ σ_full²), water vapour (∝ q), and clouds (active below the diagnosed cloud top); the column is weighted by a zenith correction factor 1 + azen (1 − μ₀)^nzen.

Fields:

• zenith_amplitude: Zenith correction amplitude (SPEEDY azen) (default NF(1))
• zenith_exponent: Zenith correction exponent (SPEEDY nzen) (default NF(2))
• absorptivity_dry_air: Absorptivity of dry air per 10⁵ Pa (default NF(0.03135))
• aerosols: Include a constant aerosol concentration (default true)
• absorptivity_aerosol: Absorptivity of aerosols per 10⁵ Pa (default NF(0.03135))
• absorptivity_water_vapor: Absorptivity of water vapour per kg/kg per 10⁵ Pa (default NF(75))
• absorptivity_cloud_base: Base cloud absorptivity per kg/kg per 10⁵ Pa (default NF(10))
• absorptivity_cloud_limit: Maximum cloud absorptivity per 10⁵ Pa (default NF(0.14))

compute_transmissivity!(
    t::AbstractVector,
    transmissivity::ConstantShortwaveTransmissivity,
    clouds,
    profile::AtmosphereProfile,
    geometry::ColumnGrid,
    surface::SurfaceState
) -> AbstractVector

Layer transmissivities under the constant-transmissivity model. Writes into t (length Nz) and returns it.

AbstractShortwaveClouds <: AbstractRadiationScheme

A diagnostic cloud model for the shortwave.

struct NoClouds <: AbstractShortwaveClouds

Trivial cloud model that reports zero cloud cover. Returned shape matches DiagnosticClouds so downstream code can be agnostic.

struct DiagnosticClouds{NF} <: AbstractShortwaveClouds

Diagnostic clouds after the Fortran SPEEDY scheme (Kucharski, Molteni & Bracco, 2006). Cloud cover is a combination of a relative-humidity term and a precipitation term; the highest layer exceeding the RH threshold sets the cloud top. An independent stratocumulus term is diagnosed at the surface from dry-static-energy stability.

Fields:

• relative_humidity_threshold_min: Relative humidity threshold for cloud cover = 0, dimensionless (default NF(0.3))
• relative_humidity_threshold_max: Relative humidity threshold for cloud cover = 1, dimensionless (default NF(1))
• specific_humidity_threshold_min: Specific humidity threshold for cloud cover in kg/kg (default NF(0.0002))
• precipitation_weight: Weight for the √precipitation term, dimensionless (default NF(0.2))
• precipitation_max: Cap on precipitation contributing to cloud cover in mm/day (default NF(10))
• cloud_albedo: Cloud albedo at CLC = 1, dimensionless (default NF(0.6))
• stratocumulus_albedo: Stratocumulus cloud albedo, dimensionless (default NF(0.5))
• stratocumulus_stability_min: Static-stability lower threshold for stratocumulus (GSES0) in J/kg (default NF(0.25))
• stratocumulus_stability_max: Static-stability upper threshold for stratocumulus (GSES1) in J/kg (default NF(0.4))
• stratocumulus_cover_max: Maximum stratocumulus cloud cover (CLSMAX), dimensionless (default NF(0.6))
• use_stratocumulus: Enable the stratocumulus parameterization (default true)
• stratocumulus_cloud_factor: Stratocumulus cloud factor (SPEEDY clfact), dimensionless (default NF(1.2))

solve_shortwave!(
    temperature_tendency::AbstractVector,
    diagnostics::ShortwaveDiagnostics{NF},
    ::NumericalRadiation.TransparentShortwave,
    profile::AtmosphereProfile,
    geometry::ColumnGrid,
    surface::SurfaceState,
    constants,
    thermodynamic;
    cloud_top_convective
)

Column transparent-atmosphere shortwave.

solve_shortwave!(
    temperature_tendency::AbstractVector,
    diagnostics::ShortwaveDiagnostics{NF},
    scheme::OneBandShortwave,
    profile::AtmosphereProfile,
    geometry::ColumnGrid,
    surface::SurfaceState,
    constants,
    thermodynamic;
    transmissivity_scratch,
    cloud_top_convective
)

Column shortwave radiative transfer for the one-band scheme.

transmissivity_scratch (length Nz) is overwritten with the layer transmissivities; pre-allocate it outside of hot loops for GPU kernels.

cloud_top_convective is the cloud top set upstream by convection or large-scale condensation; pass length(profile.temperature) + 1 if none.

rain_rate for the diagnostic cloud scheme comes from profile.rain_rate.

solve_shortwave!(
    column::RadiativeTransferColumn;
    cloud_top_convective
) -> RadiativeTransferColumn

Column shortwave radiative transfer using the scheme stored on column.

solar_declination(γ) -> Any

Solar declination [rad] from the fractional day of year γ ∈ [0, 2π] using the Spencer (1971) four-term Fourier series (accurate to ≈0.01°).

equation_of_time(γ) -> Any

Equation of time [rad] from the fractional day of year γ ∈ [0, 2π] (Spencer 1971).

cosine_solar_zenith(
    longitude,
    latitude,
    time::Dates.DateTime;
    axial_tilt,
    equinox,
    days_per_year,
    seconds_per_day
) -> Any

Cosine of the solar zenith angle at longitude and latitude [rad] and UT time. Clipped to max(μ₀, 0) so the night side returns 0.

For seasonal-only (daily-average) insolation the caller should average over a day or set time to noon and absorb the daily mean separately.

fractional_year_angle(time::Dates.DateTime) -> Float64
fractional_year_angle(
    time::Dates.DateTime,
    equinox::Dates.DateTime
) -> Float64
fractional_year_angle(
    time::Dates.DateTime,
    equinox::Dates.DateTime,
    days_per_year::Real
) -> Any

Fractional year angle γ = 2π (day − 1) / days_per_year for a DateTime.

Clausius–Clapeyron saturation specific humidity at (T, p) given ThermodynamicConstants. Returns NaN if the partial pressure is unresolvable (e.g. zero total pressure).