API Reference

adapt(
    ::Terrarium.NumberFormatAdaptor{NF},
    obj
) where NF<:Number

Adaptor that reconstructs arbitrary data structures with all numeric values converted to the specified number format NF.

Merges all of the given Variables containers into a single container.

cell_advection_timescale(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants
)
cell_advection_timescale(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology}
)
cell_advection_timescale(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology},
    timestepper::Union{Nothing, Terrarium.AbstractTimeStepper}
)

Return the minimum over the soil column of the time scales of the explicit transport terms: the advective CFL limit of the explicit first-order upwind advection of sensible heat by the soil water flux (see advective_heat_flux), and the explicit water outflow limit of the hydrology (water_outflow_timescale). The advective limit is Inf without water flow (e.g. NoFlow hydrology) and when the model's timestepper integrates internal_energy implicitly (the implicit operator is unconditionally stable). The evapotranspiration sink of an optional surface_hydrology enters the outflow limit.

cell_advection_timescale(integrator::ModelIntegrator)

Return the advective stability timescale (seconds) of the integrator, used by the Oceananigans TimeStepWizard for an advective CFL constraint. For SoilModel and LandModel this is the minimum of the timescale of the explicit upwind advection of sensible heat by the soil water flux (see advective_heat_flux), Inf without water flow and when internal_energy is integrated implicitly, and of the explicit water outflow limit (water_outflow_timescale). Models that transport nothing advectively return Inf.

cell_advection_timescale(
    state,
    grid,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Return the minimum over all cells of the advective timescale $τ = C Δz / (C_w |j_w|)$ (compute_heat_advection_timescale) for a Richards-equation hydrology, with j_w the Darcy flux through the faces of the cell.

fill_halo_regions!(field::AbstractField, state::StateVariables)

Alias for fill_halo_regions!(field, state.clock, state.inputs).

getbc(::Variable{name}, i::Integer, j::Integer, grid::AbstractGrid, clock, fields) where {name}

Implementation of Oceananigans.BoundaryConditions.getbc for variable placeholders that retrieves the input Field from fields and returns the value at the given index.

cell_diffusion_timescale(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants
)
cell_diffusion_timescale(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants,
    timestepper::Union{Nothing, Terrarium.AbstractTimeStepper}
)

Return the minimum diffusive timescale over the coupled soil processes as the minimum of the individual energy (heat conduction) and hydrology (Richards) timescales. If the model's timestepper integrates internal_energy (saturation_water_ice) implicitly (see AbstractIMEX), the energy (hydrology) timescale is Inf.

cell_diffusion_timescale(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants
)

Return the minimum thermal diffusion timescale $τ = \max(d_{\text{snow}}, d_{\text{min}})² \, C_{\text{snow}} / κ_{\text{snow}}$ over the grid for the snow-top conduction operator, mirroring the cell_diffusion_timescale pattern used for soil. Uses the same floored conduction thickness (min_snow_conduction_thickness) that enters the actual conductive-flux calculations (compute_snow_surface_heat_flux, compute_snow_basal_heat_flux), so the reported timescale reflects the stability limit of the discretization as implemented — the floor bounds it below for a vanishing pack, while a genuinely thick pack still reports its own (longer, less restrictive) timescale rather than being pinned to the floor.

cell_diffusion_timescale(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    args...
)

Return the minimum thermal diffusion timescale $τ = Δz² C / κ$ over the soil column for the SoilThermodynamics energy process, where κ is the bulk thermal conductivity (W m⁻¹ K⁻¹) and C the bulk volumetric heat capacity (J m⁻³ K⁻¹). This is the diffusive stability limit of the explicit heat conduction operator, whose effective temperature diffusivity is $α = κ / C$.

The sensible heat capacity is used (not the apparent capacity, which includes the latent heat of phase change); since the apparent capacity is never smaller, the reported timescale is conservative near the freezing front.

cell_diffusion_timescale(integrator::ModelIntegrator)

Return the minimum diffusive stability timescale $τ = Δz² / D$ (seconds) over all grid cells of the integrator, where D is the largest effective diffusivity among the model's diffusive processes. This is the diagnostic consumed by the Oceananigans TimeStepWizard to enforce a diffusive Courant–Friedrichs–Lewy (CFL) constraint of the form $Δt ≈ \mathrm{diffusive\_cfl} · τ$.

The generic fallback (for models without an implemented diffusive timescale) returns Inf, i.e. no diffusive restriction, mirroring the Oceananigans infinite_diffusion_timescale convention.

cell_diffusion_timescale(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    args...
) where NF

Return the minimum hydraulic diffusion timescale $τ = Δz² (∂θ/∂ψ) / K$ over the soil column for a Richards-equation SoilHydrology, where K is the (variably saturated) hydraulic conductivity (m s⁻¹) and ∂θ/∂ψ the specific moisture capacity (m⁻¹). This is the diffusive stability limit of the explicit Richards operator, whose effective moisture diffusivity is $D = K ∂ψ/∂θ = K / (∂θ/∂ψ)$.

As the soil saturates, $∂θ/∂ψ → 0$ and the timescale tends to zero: explicit Richards is stiff near saturation, which is precisely the regime motivating implicit timestepping. The TimeStepWizard's min_change bound keeps the step from collapsing to zero in a single adjustment.

set!(field::Field, kf::KernelFunction, fields)
set!(field::Field, kf::KernelFunction, clock, fields)

Set field by evaluating the kernel function kf at every node of field, reading from fields (and clock, for non-autonomous kernel functions). The operation is constructed on the grid of field and its location.

λnodes(grid; ...)
λnodes(grid, ::Union{Center, Face}, args...; kwargs...)

Longitudes (degrees) of the active (masked) columns of grid, as a host Vector ordered like the horizontal index of the grid. This extends the Oceananigans accessor for curvilinear grids so that λnodes/φnodes give the geographic coordinates of a horizontal position on both a LatitudeLongitudeGrid and a ColumnRingGrid; on the latter, both are indexed by column.

φnodes(grid; ...)
φnodes(grid, ::Union{Center, Face}, args...; kwargs...)

Latitudes (degrees) of the active (masked) columns of grid, as a host Vector ordered like the horizontal index of the grid. See λnodes.

serializeproperty!(file, address, grid)

Serialize a ColumnRingGrid as the RectilinearGrid it wraps. The rings and mask fields are intentionally discarded, so Fields will be read back as plain RectilinearGrids.

run!(
    integrator;
    steps,
    period,
    Δt,
    checkpointing,
    show_progress
)

Run the simulation for steps or a given time period with timestep size Δt (in seconds or Dates.Period).

Reset all Fields in state to zero.

update_state!(state::StateVariables, model::AbstractModel, inputs::InputSources; compute_tendencies = true)

Update the state for the given model and inputs; this includes calling update_inputs! and compute_auxiliary! as well as compute_tendencies!, if compute_tendencies = true.

launch!(
    grid,
    ::Union{Type{<:XY}, XY},
    kernel!,
    first_arg,
    other_args...;
    kwargs...
)

Launch kernel! over grid with the work layout given by workspec, passing grid itself as the kernel's second argument so that kernels can be written in terms of the model's grid. The work specification muist be given as a Terrarium VarDims, i.e. XY or XYZ.

Creates a TracerGas for ambient CO2 with concentration prescribed by an input variable with the given name.

ExponentialSpacing(; Δz_min, Δz_max, N)

Construct an Oceananigans ExponentialDiscretization for a soil column of N layers whose thickness grows geometrically with depth, from Δz_min at the surface to Δz_max at the bottom. The total depth of the column follows from the layer thicknesses,

D = Δz_min (ρᴺ - 1) / (ρ - 1),    ρ = (Δz_max / Δz_min)^(1 / (N - 1))

where ρ is the ratio between the thicknesses of successive layers. Oceananigans parameterizes the same geometric family by the e-folding scale of its exponential coordinate mapping, which for the spacings above is D / (N log ρ). Note that Oceananigans orders cell interfaces bottom-up, so the first cell of the returned discretization is the thick bottom layer.

For a column of constant layer thickness, use UniformSpacing instead; the exponential mapping is degenerate when Δz_min == Δz_max.

using Terrarium

z = ExponentialSpacing(Δz_min = 0.05, Δz_max = 100.0, N = 50)
(num_layers(z), round(z.faces[1], digits = 3), z.faces[end])

# output
(50, -695.654, 0.0)

Alias for PrescribedGradient representing a Neumann-type zero pressure gradient at the bottom of the soil column, thereby allowing free drainage of water.

Alias for FluxBoundaryCondition on internal_energy with name geothermal_heat_flux representing the geothermal heat flux at the bottom boundary of the soil column.

HomogeneousSoilStratigraphy(
    ::Type{NF};
    texture::SoilTexture,
    porosity::Terrarium.AbstractSoilPorosity
) where NF

Convenience constructor that creates a SoilStratigraphy with a single ConstantSoilHorizon.

Alias for NoFlux representing a zero-flux bottom boundary condition for water flow (prognostic variable saturation_water_ice).

Alias for PrescribedFlux with name infiltration representing liquid water infiltration at the soil surface.

Alias for a discrete-form FluxBoundaryCondition on saturation_water_ice that normalizes the physical infiltration flux (m/s of water depth) by the top-of-soil porosity at evaluation time, since saturation_water_ice is the dimensionless saturation (VWC / porosity). See saturation_infiltration_bc.

Alias for FluxBoundaryCondition on density_soc with name litter representing litter fall inputs to top layer soil organic carbon.

Alias for ValueBoundaryCondition on top temperature (in °C) with the given variable name.

Alias for ValueBoundaryCondition on top temperature (in °C) with the given variable name.

Creates a Callback for ProgressReporter with schedule IterationInterval(1).

SoilGridsStratigraphy(
    ::Type{NF};
    porosity::Terrarium.AbstractSoilPorosity,
    horizon1::Terrarium.AbstractSoilHorizon,
    horizon2::Terrarium.AbstractSoilHorizon,
    horizon3::Terrarium.AbstractSoilHorizon,
    horizon4::Terrarium.AbstractSoilHorizon,
    horizon5::Terrarium.AbstractSoilHorizon,
    horizon6::Terrarium.AbstractSoilHorizon
) where NF -> SoilStratigraphy{NF, 6, <:NTuple{6, Terrarium.AbstractSoilHorizon{NF}}} where NF

Convenience constructor that creates a SoilStratigraphy with six horizons corresponding to the six depth intervals defined by the SoilGridsV2 dataset. The porosity parameterization is applied to all six default horizons; pass explicit horizon* arguments to override individual horizons (in which case porosity does not apply to those).

Alias for FluxBoundaryCondition on internal_energy with name soil_heat_flux representing the net heat flux into the top of the soil column. Without snow this equals the surface energy balance ground_heat_flux; with snow it is the (blended) conductive flux across the snow base.

Creates a NamedTuple from the given tracer gas types.

UniformSpacing(; Δz, N)

Construct the cell interfaces of a soil column of N layers of constant thickness Δz, spanning a total depth of N * Δz with the surface at z = 0. The number format of the returned interfaces follows that of Δz, so it also determines the default number format of a grid built from them.

This is a convenience for the uniformly spaced case, which Oceananigans expresses as a plain range of cell interfaces; UniformSpacing(Δz = 0.1, N = 10) is equivalent to range(-10 * 0.1, 0, length = 11). See ExponentialSpacing for a column whose layers grow with depth.

using Terrarium

z = UniformSpacing(Δz = 0.1, N = 10)
(num_layers(z), first(z), last(z))

# output
(10, -1.0, 0.0)

Generation-time helper for get_fields. Given the type Vars of a tuple of AbstractVariables and Namespaces, builds an expression that retrieves all matching fields from state_ex (an expression evaluating to the state container) and vars_ex (an expression evaluating to the variable tuple). The recursion over namespaces is performed here, at expansion time, so that the generated body for get_fields contains no self-call. This is what makes the method type stable: a runtime self-recursive get_fields would otherwise trigger inference's recursion limiting and widen the return type of the nested namespace lookups to an abstract NamedTuple.

add_scaled!(out, field, factor)

Add factor times field to out over the interior of their common grid, out ← out + factor * field, e.g. to form or remove an explicit increment. Fields without a vertical location (e.g. snow variables) are updated once per column. Not exported.

adjust_saturation_profile!(state, grid, hydrology::SoilHydrology, runoff = nothing, strat = nothing, bgc = nothing, energy = nothing)

Redistribute over- and undersaturation in the soil columns of state (see redistribute_saturation_profile!). With runoff = nothing, excess water reaching the surface is discarded; otherwise it is routed into the surface_excess_water pool of runoff as a water depth (m), using the porosity of the topmost layer from strat and bgc.

With an optional SoilThermodynamics energy process, the sensible heat of the water moves with it (compute_redistributed_heat!), so the heat of the excess water leaves the column; this requires strat and bgc.

adjust_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    ::Nothing,
    energy,
    strat,
    bgc
) where NF

Kernel function that adjusts the saturation profile at i, j and discards any excess water reaching the surface. Used for standalone soil hydrology, where no surface runoff process owns a water pool.

adjust_saturation_profile!(
    out,
    i,
    j,
    grid,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    runoff::Terrarium.AbstractSurfaceRunoff,
    energy,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Kernel function that adjusts the saturation profile at i, j and routes any excess water reaching the surface into the surface_excess_water pool owned by the given surface runoff process, as a water depth (m): the excess saturation times the thickness and porosity (porosity_top) of the topmost layer.

advective_heat_flux(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Advective heat flux (W/m², positive upward) at the cell face i, j, k. The default is zero, i.e. no heat is transported by a moving medium; see the method for soil water flow following the Richards equation.

advective_heat_flux(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    args...
) where NF

Upwind advective heat flux C_w T_up j_w (W/m², positive upward) at face i, j, k, the first-order upwind discretization of C_w T j_w, with j_w (m/s) the water_flux of the step (compute_advective_water_flux) and C_w the water heat capacity of the SoilThermalProperties, consistent with the energy closure. T_up is the temperature of cell k for j_w < 0 and of cell k - 1 otherwise, clamped to the column at the boundary faces. Only the sensible heat of liquid water is advected; ice is immobile.

Under Implicit time stepping, only the boundary inflow is explicit (compute_boundary_advective_energy_tendency); the rest is part of the implicit heat operator (compute_implicit_advective_velocity!). The heat of water entering through flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) is added by compute_boundary_advected_heat!.

aerodynamic_resistance(i, j, grid, fields, atmos::PrescribedAtmosphere)

Compute the aerodynamic resistance (inverse conductance) at grid cell i, j.

aerodynamic_resistance(
    i,
    j,
    grid,
    fields,
    atmos::Terrarium.AbstractAtmosphere,
    evapotranspiration::PALADYNCanopyEvapotranspiration{NF},
    ::Terrarium.AbstractVegetation
) where NF

Compute the aerodynamic resistance between the ground and canopy as a function of LAI and SAI.

air_density(i, j, grid, fields, atmos::AbstractAtmosphere, constants::PhysicalConstants)

Compute density (kg m⁻³) of a parcel of air under current atmospheric conditions.

air_pressure(i, j, grid, fields, ::PrescribedAtmosphere)

Retrieve or compute the air pressure at the current time step.

air_temperature(i, j, grid, fields, ::PrescribedAtmosphere)

Retrieve or compute the air temperature at the current time step.

albedo(i, j, grid, fields, albedo::ConstantAlbedo)

Return the surface albedo at grid point i, j.

albedo(i, j, grid, fields, ::AbstractAlbedo)

Return the current albedo at the given indices.

ambient_co2(i, j, grid, fields, ::PrescribedAtmosphere)

Return the current prescribed ambient CO2 concentration level.

auxiliary(name, loc; ...)
auxiliary(name, loc, ctor; ...)
auxiliary(
    name,
    loc,
    ctor,
    params;
    units,
    eltype,
    closure,
    bounds,
    desc
)

Convenience constructor method for AuxiliaryVariable.

auxiliary_fields(state, components...)

Retrieves all Fields from state corresponding to auxiliary variables defined on the given components.

Helper method that selects only auxiliary variables declared on obj.

canopy_domain(::Oceananigans.Grids.AbstractGrid)

Return the spatial discretization of the canopy domain of grid, or nothing if the canopy is not vertically resolved. Grids which are not land grids never resolve a canopy domain.

canopy_evaporation_conductance(
    ::PALADYNCanopyEvapotranspiration,
    f_can,
    rₐ
)

Compute the canopy evaporation vapor conductance [m/s] from the current canopy saturation fraction f_can and aerodynamic resistance rₐ.

canopy_water(i, j, grid, fields, ::AbstractCanopyInterception)

Compute or retrieve the current canopy water storage [m].

celsius_to_kelvin(c::ThermodynamicConstants, T)

Convert the given temperature in °C to Kelvin based on the constant temperature_water_freeze.

check_domain_compatibility(ground, domain, name)

Check that the domain grid named name is compatible with the ground domain grid, i.e. that both share the same horizontal discretization, numeric type, and architecture. Throws an ArgumentError if not. Note that this is a host-side check which must never be invoked from a kernel.

Check for variables/namespaces with duplicate names and raise an error if duplicates are detected. Not type stable.

checkfinite!(field)
checkfinite!(field, name)

Check whether the given field has any NaN or Inf values and raise an error if NaNs are detected.

closure!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants
)
closure!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology}
)

Compute the forward closure mapping for soil hydrology and energy, in that order.

An optional surface_hydrology process may be supplied (by the coupled LandModel) so that excess water removed from an oversaturated soil surface is routed into the surface_excess_water pool owned by its runoff scheme (get_surface_runoff). Without it (standalone soil), the excess is discarded.

closure!(
    state,
    grid,
    closure::SnowEnergyTemperatureClosure,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    args...
)

Forward closure: recover snow_temperature and snow_liquid_fraction from the prognostic snow_energy.

closure!(state, grid, closure::AbstractClosureRelation, process, args...)

Apply closure for process on the given grid with additional implementation-specific args.

closure!(state, grid, proc::AbstractProcess, args...)

Apply the forward closure mappings for the process proc on the given grid with additional implementation-specific args defined by the coupling interface for the process type. The default implementation calls invclosure! with args for each closure returned by closures.

closure!(state, model::AbstractModel)

Apply all closure relations defined for the given model, i.e. update the variables that are determined by the prognostic state. This is also where a model enforces constraints on its prognostic state that are not part of the tendencies, such as the redistribution of oversaturation in the soil (adjust_saturation_profile!) or the non-negativity of the snow water equivalent (enforce_snow_constraints!). Timesteppers call it at the end of each step (and between the stages of multi-stage schemes), after all prognostic variables have been updated.

closure!(
    state,
    grid,
    closure::SoilSaturationPressureClosure,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil
) where NF
closure!(
    state,
    grid,
    closure::SoilSaturationPressureClosure,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    args...
) where NF

Computes pressure_head $Ψ = ψm + ψz + ψh$ from the current saturation_water_ice state, after redistributing over- and undersaturation (adjust_saturation_profile!) and updating the water table.

An optional runoff process routes excess water removed from an oversaturated soil surface into the runoff-owned surface_excess_water pool; without it the excess is discarded. With SoilThermodynamics, the sensible heat of redistributed water is moved in internal_energy and that of the excess is removed; the energy closure, which runs after this one in closure! of SoilEnergyWaterCarbon, then recomputes the temperature.

closure_fields(state, components...)

Retrieves all Fields from state corresponding to closure variables defined on the given components.

Helper method that selects only closure (auxiliary) variables declared on obj.

closures(proc::AbstractProcess)

Return a tuple of AbstractClosureRelations defined by the given processes type. Note that this is a type-stable, @generated function that is compiled for each argument type.

compute_APAR(
    photo::LUEPhotosynthesis{NF},
    traits::PlantTraits{NF},
    swdown::NF,
    LAI::NF
) where NF

Computes absorbed PAR limited by the fraction of PAR assimilated at ecosystem level APAR [mol/m²/s], [21, Eq. (62)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_Ag(
    photo::LUEPhotosynthesis{NF},
    c_1::NF,
    c_2::NF,
    APAR::NF,
    Vc_max::NF,
    β::NF
) where NF

Computes the gross photosynthesis rate Ag [kgC/m²/s], [22, Eq. (2)].

References

• [22] Haxeltine & Prentice, Global Biogeochemical Cycles (1996)

compute_C_veg_tend(vegcarbon_dynamics, traits, LAI_b, NPP)

Computes the C_veg tendency based on NPP and the balanced Leaf Area Index LAI_b, [21, Eq. (72)]

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_GPP(
    ::LUEPhotosynthesis{NF},
    An::NF,
    Rd::NF
) where NF

Compute Gross Primary Production from net assimilation An and leaf respiration Rd [kgC/m²/s].

compute_JE_JC(
    photo::LUEPhotosynthesis{NF},
    c_1::NF,
    c_2::NF,
    APAR::NF,
    Vc_max::NF
) where NF -> Tuple{Any, Any}

Computes the PAR-limited and the rubisco-activity-limited photosynthesis rates JE and JC [kgC/m²/s], [22, Eqs. (3, 5)].

References

• [22] Haxeltine & Prentice, Global Biogeochemical Cycles (1996)

compute_NPP(autoresp, GPP, Ra)

Computes Net Primary Productivity NPP as the difference between Gross Primary Production GPP and autotrophic respiration Ra in [kgC/m²/s].

compute_PAR(
    photo::LUEPhotosynthesis{NF},
    swdown::NF
) where NF

Computes NET Photosynthetically Active Radiation PAR [mol/m²/s].

compute_Ra(
    autoresp,
    vegcarbon_dynamics,
    traits,
    T_air,
    T_soil,
    Rd,
    phen,
    C_veg,
    GPP
)

Computes autotrophic respiration Ra as the sum of maintenance respiration Rm and growth respiration Rg in [kgC/m²/s].

compute_Rd(
    photo::LUEPhotosynthesis,
    Vc_max::NF,
    β::NF
) where NF

Computes the leaf respiration rate Rd [kgC/m²/s], [22, Eq. (10)] and [21, Eq. (10)].

References

• [22] Haxeltine & Prentice, Global Biogeochemical Cycles (1996)
• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_Rg(autoresp, GPP, Rm)

Computes growth respiration Rg in [kgC/m²/s].

compute_Rm(
    autoresp,
    vegcarbon_dynamics,
    traits,
    T_air,
    T_soil,
    Rd,
    phen,
    C_veg
)

Computes maintenance respiration Rm in [kgC/m²/s].

compute_Vc_max(
    photo::LUEPhotosynthesis{NF},
    c_1::NF,
    PAR::NF,
    Kc::NF,
    Ko::NF,
    Γ_star::NF,
    pres_i::NF,
    pres_O2::NF
) where NF

Computes the maximum carboxylation capacity Vc_max [kgC/m²/s], following the coordination hypothesis (acclimation), see [23] Box 2. Note: this is not the same formula as in the [21] paper, this implementation is taken from the PALADYN code

References

• [23] Harrison et al., New Phytologist (2021)
• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_advective_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Compute the part -∂z(q_adv) of the internal energy tendency at i, j, k that is due to the divergence of the advective_heat_flux q_adv, through all faces of the cell. When the heat fluxes are integrated implicitly, only the part due to the boundary faces is treated explicitly (compute_boundary_advective_energy_tendency).

compute_advective_water_flux(i, j, k, grid, fields)

Vertical water flux j_w (m/s, positive upward) at face i, j, k that advects heat: the water_flux of the step, i.e. the Darcy flux at tⁿ (boundary faces, and interior faces under explicit hydrology) or of the implicit update (see compute_water_flux and compute_implicit_water_flux).

compute_albedo!(
    out,
    i,
    j,
    grid,
    fields,
    albedo::DiagnosticAlbedo,
    vegetation::Union{Nothing, Terrarium.AbstractVegetation},
    snow::Union{Nothing, Terrarium.AbstractSnow},
    args...
)

compute_albedo(
    i,
    j,
    grid,
    fields,
    snow::Terrarium.AbstractSnow
)

compute_albedo(
    i,
    j,
    grid,
    fields,
    albedo::Terrarium.ConstantSnowAlbedo
)

compute_assimilation_factors(
    photo::LUEPhotosynthesis{NF},
    constants::MaterialConstants{NF},
    Γ_star::NF,
    T_stress::NF,
    Kc::NF,
    Ko::NF,
    pres_i::NF,
    pres_O2::NF
) where NF -> Tuple{Any, Any}

Computes factors for light-limited c_1 [kgC/mol] and RuBisCO-limited c_2 [dimensionless] assimilation. Follows [21, Eqs. (C4-C5)] and [22].

• c1: quantum efficiency × temperature factor × carbon mass / denominator, used in JE = c1 × APAR
• c2: dimensionless coefficient relating enzyme capacity to assimilation, used in JC = c2 × Vc_max

References

• [22] Haxeltine & Prentice, Global Biogeochemical Cycles (1996)
• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)H

compute_autotrophic_respiration(
    i, j, grid, fields,
    autoresp::AbstractAutotrophicRespiration,
    vegcarbon::AbstractVegetationCarbonDynamics,
    atmos::AbstractAtmosphere,
    args...
)

Cell-level autotrophic respiration computation. Implementations should compute autotrophic respiration and related diagnostics (e.g. NPP) for the given cell.

compute_autotrophic_respiration!(
    out,
    i,
    j,
    grid,
    fields,
    autoresp::Terrarium.AbstractAutotrophicRespiration,
    args...
)

Mutating wrapper for compute_autotrophic_respiration that stores the results in out.

compute_autotrophic_respiration(
    i,
    j,
    grid,
    fields,
    autoresp::PALADYNAutotrophicRespiration,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    phenology::Terrarium.AbstractPhenology,
    traits::PlantTraits,
    atmos::Terrarium.AbstractAtmosphere
) -> Tuple{Any, Any}

Compute autotrophic respiration following the scheme of [21].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_auxiliary!(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants
)
compute_auxiliary!(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    timestepper::Union{Nothing, Terrarium.AbstractTimeStepper}
)

Diagnose the snow geometric properties (snow_depth, snow_cover_fraction) from the current snow water equivalent and bulk density, and the meltwater_flux from the liquid fraction if snow_water_equivalent is integrated explicitly by timestepper. Under implicit integration, meltwater_flux holds the drainage of the last implicit step instead (see implicit_update!).

compute_auxiliary!(
    state,
    grid,
    evaporation::BareGroundEvaporation,
    ::NoCanopyInterception,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_auxiliary!(
    state,
    grid,
    evaporation::BareGroundEvaporation,
    ::NoCanopyInterception,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil}
)
compute_auxiliary!(
    state,
    grid,
    evaporation::BareGroundEvaporation,
    ::NoCanopyInterception,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil},
    snow::Union{Nothing, Terrarium.AbstractSnow}
)

compute_auxiliary!(state, model::AbstractModel)

Compute updates to all auxiliary variables based on the current prognostic state of the model.

compute_auxiliary!(state, grid, process::AbstractProcess, args...)

Compute all auxiliary state variables for the given process on grid. Implementations of AbstractProcess may define additional args that correspond to different process coupling interfaces.

compute_auxiliary!(
    state,
    grid,
    runoff::DirectSurfaceRunoff,
    canopy_interception::Terrarium.AbstractCanopyInterception,
    soil::Terrarium.AbstractSoil
)
compute_auxiliary!(
    state,
    grid,
    runoff::DirectSurfaceRunoff,
    canopy_interception::Terrarium.AbstractCanopyInterception,
    soil::Terrarium.AbstractSoil,
    snow::Union{Nothing, Terrarium.AbstractSnow},
    args...
)

compute_auxiliary!(
    state,
    grid,
    veg::PrognosticVegetation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_auxiliary!(
    state,
    grid,
    veg::PrognosticVegetation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil},
    args...
)

Compute auxiliary variables for all vegetation component processes based on the given atmospheric inputs defined by atmos and (optionally) soil state. If soil = nothing, stress factors due to soil temperature and moisture availability will be ignored.

compute_auxiliary!(state, grid, albedo::DiagnosticAlbedo)
compute_auxiliary!(
    state,
    grid,
    albedo::DiagnosticAlbedo,
    vegetation::Union{Nothing, Terrarium.AbstractVegetation}
)
compute_auxiliary!(
    state,
    grid,
    albedo::DiagnosticAlbedo,
    vegetation::Union{Nothing, Terrarium.AbstractVegetation},
    snow::Union{Nothing, Terrarium.AbstractSnow},
    args...
)

Diagnose the blended surface albedo and emissivity, weighting the bare-ground, vegetation, and snow values by the vegetated area fraction of the optional vegetation component and the snow-covered area fraction of the optional snow component.

compute_auxiliary!(
    state,
    grid,
    hydrology::SurfaceHydrology,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_auxiliary!(
    state,
    grid,
    hydrology::SurfaceHydrology,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil}
)
compute_auxiliary!(
    state,
    grid,
    hydrology::SurfaceHydrology,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil},
    vegetation::Union{Nothing, Terrarium.AbstractVegetation}
)
compute_auxiliary!(
    state,
    grid,
    hydrology::SurfaceHydrology,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil},
    vegetation::Union{Nothing, Terrarium.AbstractVegetation},
    snow::Union{Nothing, Terrarium.AbstractSnow},
    args...
)

compute_auxiliary!(
    state,
    grid,
    veg::PrescribedVegetation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_auxiliary!(
    state,
    grid,
    veg::PrescribedVegetation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil},
    args...
)

Compute auxiliary variables for all vegetation component processes based on the given atmospheric inputs defined by atmos and (optionally) soil state. If soil = nothing, stress factors due to soil temperature and moisture availability will be ignored.

compute_auxiliary!(state, grid, seb::SurfaceEnergyBalance)
compute_auxiliary!(
    state,
    grid,
    seb::SurfaceEnergyBalance,
    vegetation::Union{Nothing, Terrarium.AbstractVegetation}
)
compute_auxiliary!(
    state,
    grid,
    seb::SurfaceEnergyBalance,
    vegetation::Union{Nothing, Terrarium.AbstractVegetation},
    snow::Union{Nothing, Terrarium.AbstractSnow},
    args...
)

compute_auxiliary!(
    state,
    grid,
    rad::DiagnosedRadiativeFluxes,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    consts::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

compute_auxiliary!(
    state,
    grid,
    tur::DiagnosedTurbulentFluxes,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

compute_auxiliary!(
    state,
    grid,
    paw::FieldCapacityLimitedPAW,
    soil::Terrarium.AbstractSoil,
    args...
)

compute_auxiliary!(
    state,
    grid,
    photo::LUEPhotosynthesis,
    stomcond::Terrarium.AbstractStomatalConductance,
    traits::PlantTraits,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

compute_auxiliary!(
    state,
    grid,
    stomcond::MedlynStomatalConductance,
    traits::PlantTraits,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

compute_auxiliary!(
    state,
    grid,
    autoresp::PALADYNAutotrophicRespiration,
    vegcarbon::Terrarium.AbstractVegetationCarbonDynamics,
    phenology::Terrarium.AbstractPhenology,
    traits::PlantTraits,
    atmos::Terrarium.AbstractAtmosphere
)

compute_auxiliary!(
    state,
    grid,
    evapotranspiration::PALADYNCanopyEvapotranspiration,
    interception::Terrarium.AbstractCanopyInterception,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Terrarium.AbstractSoil,
    vegetation::Terrarium.AbstractVegetation,
    args...
)

compute_auxiliary!(
    state,
    grid,
    canopy_interception::PALADYNCanopyInterception,
    atmos::Terrarium.AbstractAtmosphere
)

compute_auxiliary!(
    state,
    grid,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits,
    args...
)

compute_auxiliary!(
    state,
    grid,
    phenol::PALADYNPhenology,
    vegcarbon::PALADYNCarbonDynamics,
    atmos::Terrarium.AbstractAtmosphere
)

compute_auxiliary!(
    state,
    grid,
    veg_dynamics::PALADYNVegetationDynamics,
    args...
)

compute_auxiliary!(
    state,
    grid,
    rad::PrescribedRadiativeFluxes,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

compute_auxiliary!(
    state,
    grid,
    soc::SinglePoolSoilCarbon,
    soil::Terrarium.AbstractSoil,
    args...
)

Compute the soil organic carbon respiration rate in each soil layer.

compute_auxiliary!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants
)

Compute auxiliary variables for soil energy, water, and carbon state variables on grid based on the given values in constants.

compute_auxiliary!(
    state,
    grid,
    hydrology::SoilHydrology,
    soil::Terrarium.AbstractSoil,
    args...
)

compute_auxiliary!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    args...
)

compute_auxiliary!(
    state,
    grid,
    ghf::Terrarium.AbstractGroundHeatFlux,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    args...
)

compute_auxiliary!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    args...
) where NF

compute_balanced_leaf_area_index(
    vegcarbon_dynamics,
    traits,
    C_veg
)

Computes LAI_b, the balanced Leaf Area Index based on the vegetation carbon pool C_veg (assuming with bwl = 1), [21, Eqs. (76-79)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_boundary_advected_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    saturation_tendency,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Internal energy tendency C_w T φ G_s (W/m³) at i, j, k of water entering (G_s > 0) or leaving (G_s < 0) the cell through flux boundary conditions on saturation_water_ice, with G_s = saturation_tendency[i, j, k], porosity φ, and cell temperature T. Infiltration thus enters at the temperature of the top cell.

compute_boundary_advected_heat!(
    state,
    grid,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Add the sensible heat of water crossing the column boundaries through flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) to the internal energy tendency (compute_boundary_advected_energy_tendency). Must be called right after the hydrology boundary conditions, while the saturation tendency holds only their contributions. No-op without water flow.

compute_boundary_advective_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Compute the part of the internal energy tendency ∂U∂t at i, j, k due to the advective heat flux through the boundary faces of the column that remains explicit when the interior fluxes are integrated implicitly (see BackwardEuler). The default is zero, i.e. no heat is transported by a moving medium; see the method for soil water flow following the Richards equation.

compute_boundary_advective_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    args...
) where NF

Explicit advective energy tendency at i, j, k under Implicit time stepping: the heat of water entering the column through the boundary faces with the Darcy flux at tⁿ, at the temperature of the boundary cell. Outflow is implicit (compute_implicit_advection_tendency), so that it cannot remove more heat than a cell holds.

compute_boundary_conditions!(state, model::AbstractModel)

Compute all internal and external boundary conditions for spatially explicit prognostic and auxiliary variables defined on the model. This should typically be invoked after compute_auxiliary! but before compute_tendencies!.

compute_boundary_conditions!(state, grid, process::AbstractProcess, args...)

Compute boundary conditions for all spatially explicit prognostic and auxiliary variables defined by the given process on grid.

compute_boundary_conditions!(
    state,
    grid,
    ::SinglePoolSoilCarbon
)

Fill the halo regions of density_soc and apply its boundary conditions (e.g. LitterfallFlux) to the soil organic carbon tendency.

compute_boundary_conditions!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon
)

Compute boundary conditions (and halo regions) for soil energy and hydrology. The sensible heat of water entering or leaving the column through flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) is added to the energy tendency by compute_boundary_advected_heat!, after the hydrology boundary conditions (the saturation tendency then holds only their contribution).

compute_boundary_conditions!(
    state,
    grid,
    ::SoilThermodynamics
)

compute_boundary_conditions!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Fill the pressure_head halo and apply the top/bottom flux boundary conditions for saturation_water_ice. The fields passed to those boundary conditions are built explicitly here (rather than using state.inputs, which only carries declared input variables) so that a discrete-form condition like saturation_infiltration_bc can read auxiliary fields (infiltration) and per-horizon namespaces (via strat/bgc) that only exist once the full state does. infiltration is only merged in when present, since standalone (non-LandModel) configurations may not declare it at all.

compute_boundary_darcy_flux_divergence(
    i,
    j,
    k,
    grid,
    fields
)

Part of the Darcy flux divergence at i, j, k due to the two boundary faces of the column (bottom face of cell k = 1, top face of cell k = Nz), omitting all interior faces:

(δ_{k,Nz} j_{w,k+½} - δ_{k,1} j_{w,k-½}) / Δz_k,

with j_w the darcy_flux. This is the explicit part when the interior fluxes are integrated implicitly (see BackwardEuler). The boundary fluxes are computed from the pressure_head halo: they vanish for the default (no-flux) conditions and carry the prescribed flux for a gradient condition on pressure_head. Flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) are added separately by compute_boundary_conditions!.

compute_boundary_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Internal energy tendency ∂U∂t at i, j, k from the diffusive_heat_flux q through the boundary faces of the column only (bottom of k = 1, top of k = Nz):

∂U∂t = -(δ_{k,Nz} q_{k+½} - δ_{k,1} q_{k-½}) / Δz_k.

Zero for flux and no-flux conditions, whose halos have zero gradient (flux conditions are added separately); nonzero, e.g., for a gradient condition on temperature.

compute_canopy_interception(
    canopy_interception::PALADYNCanopyInterception{NF},
    precip,
    LAI,
    SAI
) where NF

Compute I_can, the canopy rain interception, following [21, Eq. (42)].

References

• [21] Willeit and Ganopolski, Geoscientific Model Development (2016)

compute_canopy_saturation_fraction(
    canopy_interception::PALADYNCanopyInterception{NF},
    W_can,
    LAI,
    SAI
) where NF

Compute the canopy saturation fraction as W_can / W_can_max.

compute_canopy_water_removal(
    canopy_interception::PALADYNCanopyInterception{NF},
    W_can
) where NF

Compute the canopy water removal rate as W_can / τw.

compute_canopy_water_tendency(
    ::PALADYNCanopyInterception{NF},
    I_can,
    E_can,
    R_can
) where NF

Compute the W_can tendency and removal rate following [21, Eq. (41)].

References

• [21] Willeit and Ganopolski, Geoscientific Model Development (2016)

compute_constant_saturation(
    i,
    j,
    k,
    grid,
    fields,
    ::ConstantSaturation
)

Initial saturation at cell (i, j, k) taken from the initial_saturation field.

compute_constant_temperature(
    i,
    j,
    k,
    grid,
    fields,
    ::ConstantSoilTemperature
)

Initial temperature at cell (i, j, k) taken from the initial_surface_temperature (°C) field.

compute_darcy_flux_divergence(
    i,
    j,
    k,
    grid,
    fields,
    ::Explicit
)

Divergence ∂j_w/∂z (1/s) of the darcy_flux j_w at cell i, j, k: over all faces under Explicit time stepping, and over the boundary faces of the column only (see compute_boundary_darcy_flux_divergence) under Implicit.

compute_emissivity(
    i,
    j,
    grid,
    fields,
    snow::Terrarium.AbstractSnow
)

compute_emissivity(
    i,
    j,
    grid,
    fields,
    albedo::Terrarium.ConstantSnowAlbedo
)

compute_energy_tendencies!(
    tendencies,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    ::Explicit,
    args...
)

Add the explicitly integrated part of the internal energy tendency at i, j, k when internal_energy is integrated explicitly, i.e. the full tendency (compute_energy_tendency).

compute_energy_tendencies!(
    tendencies,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    ::Implicit,
    args...
)

Add the explicitly integrated part of the internal energy tendency at i, j, k when internal_energy is integrated implicitly: the interior fluxes belong to the implicit operator, so only the diffusive fluxes through the boundary faces (compute_boundary_energy_tendency) and the heat of water entering through them (compute_boundary_advective_energy_tendency) remain.

compute_energy_tendency(i, j, k, grid, ::AbstractThermodynamics, args...)

Compute the internal energy tendency ∂U∂t at index i, j, k.

compute_energy_tendency(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Compute the internal energy tendency ∂U∂t as the divergence of the diffusive and advective heat fluxes,

∂U∂t = -∂z(q_diff + q_adv),

where q_diff is the diffusive_heat_flux and q_adv the advective_heat_flux (zero unless the soil water flows).

compute_evapotranspiration_conductances(
    i,
    j,
    grid,
    fields,
    evaporation::BareGroundEvaporation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_evapotranspiration_conductances(
    i,
    j,
    grid,
    fields,
    evaporation::BareGroundEvaporation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil}
)

Compute the skin-driven ground evaporation vapor conductance β/rₐ at grid cell i, j for the given bare-ground evaporation scheme. This conductance is independent of the skin temperature.

compute_evapotranspiration_conductances!(
    out,
    i,
    j,
    grid,
    fields,
    evaporation::BareGroundEvaporation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_evapotranspiration_conductances!(
    out,
    i,
    j,
    grid,
    fields,
    evaporation::BareGroundEvaporation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Union{Nothing, Terrarium.AbstractSoil}
)

Compute and store the skin-driven ground evaporation vapor conductance on grid for the given bare-ground evaporation scheme.

compute_evapotranspiration_conductances!(
    out,
    i,
    j,
    grid,
    fields,
    evapotranspiration::PALADYNCanopyEvapotranspiration,
    interception::Terrarium.AbstractCanopyInterception,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Terrarium.AbstractSoil,
    vegetation::Terrarium.AbstractVegetation,
    args...
)

Compute and store the skin-driven vapor conductances on grid for the given scheme evapotranspiration and process dependencies.

compute_evapotranspiration_conductances(
    i,
    j,
    grid,
    fields,
    evapotranspiration::PALADYNCanopyEvapotranspiration,
    interception::Terrarium.AbstractCanopyInterception,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    soil::Terrarium.AbstractSoil,
    vegetation::Terrarium.AbstractVegetation,
    args...
) -> Tuple{Any, Any, Any}

Compute the skin-driven vapor conductances (ground_evaporation_conductance, canopy_evaporation_conductance, and transpiration_conductance) at grid cell i, j for the given scheme evapotranspiration and process dependencies. These conductances are independent of the skin temperature.

Compute transpiration, evaporation_ground, and evaporation_canopy fluxes on grid for the given scheme evapotranspiration and process dependencies.

compute_f_temp(autoresp, T_air, T_soil)

Computes temperature factors f_temp_air and f_temp_soil for autotrophic respiration.

compute_gdd_tendency!(
    tend,
    i,
    j,
    grid,
    fields,
    phenol::PALADYNPhenology,
    atmos::Terrarium.AbstractAtmosphere
)

Mutating wrapper for compute_gdd_tendency that stores the growing-degree-day tendency in tend.

compute_gdd_tendency(phenol, gdd, T_air)

Compute the growing-degree-day tendency d(gdd)/dt [K⋅day/s] at a single grid point.

Heat above the base temperature is accumulated (in degree-days per second) and, below the base temperature, the accumulator relaxes toward zero with timescale gdd_relaxation_time. The relaxation provides a continuous surrogate for the discrete annual reset of the original scheme.

compute_ground_heat_flux!(
    out,
    i,
    j,
    grid,
    fields,
    ghf::DiagnosedGroundHeatFlux,
    skinT::Terrarium.AbstractSkinTemperature,
    seb::Terrarium.AbstractSurfaceEnergyBalance
)

Per-cell mutating variant used by the fused surface-energy-balance kernel: store the ground heat flux into the auxiliary output field out.

compute_ground_heat_flux!(
    out,
    i,
    j,
    grid,
    fields,
    ::PrescribedGroundHeatFlux,
    ::Terrarium.AbstractSkinTemperature,
    ::Terrarium.AbstractSurfaceEnergyBalance
)

The prescribed ground heat flux is an input field, so the fused surface-energy-balance kernel leaves it untouched.

compute_ground_heat_flux!(
    state,
    grid,
    ::PrescribedGroundHeatFlux,
    args...
)

The ground heat flux is a prescribed input variable, so there is nothing to diagnose.

compute_ground_heat_flux!(
    state,
    grid,
    ghf::Terrarium.AbstractGroundHeatFlux,
    skinT::Terrarium.AbstractSkinTemperature,
    seb::Terrarium.AbstractSurfaceEnergyBalance
)

Compute and store ground_heat_flux on grid, dispatching on ghf and skinT to the type-specific kernel function. For PrescribedGroundHeatFlux this is a no-op, since the field is an input.

compute_ground_heat_flux(
    i,
    j,
    grid,
    fields,
    ::DiagnosedGroundHeatFlux,
    skinT::ImplicitSkinTemperature,
    ::Terrarium.AbstractSurfaceEnergyBalance
)

Compute the conductive ground heat flux from the current skin_temperature and ground_temperature.

compute_ground_heat_flux(
    i,
    j,
    grid,
    fields,
    ghf::DiagnosedGroundHeatFlux,
    ::PrescribedSkinTemperature,
    seb::Terrarium.AbstractSurfaceEnergyBalance
)

For PrescribedSkinTemperature, set the ground heat flux directly to the demand, i.e. G₀ = R_net + H_s + H_l.

compute_ground_heat_flux_demand(ghf::AbstractGroundHeatFlux, args...)

Return the ground heat flux demanded by the rest of the surface energy balance, which the skin temperature schemes invert their conduction relations against.

compute_ground_heat_flux_demand(
    i,
    j,
    grid,
    fields,
    ghf::DiagnosedGroundHeatFlux,
    seb::Terrarium.AbstractSurfaceEnergyBalance
)

Compute the ground heat flux demand at grid cell i, j: the flux implied by the radiative budget and the turbulent fluxes, G = R_net + H_s + H_l. This is the quantity the implicit skin temperature solve inverts its conduction relation against. The radiative and turbulent flux terms are read through the accessors of the seb sub-processes that own them.

compute_ground_heat_flux_demand(
    i,
    j,
    grid,
    fields,
    ghf::PrescribedGroundHeatFlux,
    ::Terrarium.AbstractSurfaceEnergyBalance
)

When the ground heat flux is prescribed, the demand is the prescribed flux; the surface energy budget is not re-closed from the radiative and turbulent terms.

compute_ground_heat_flux_demand(
    ::DiagnosedGroundHeatFlux,
    R_net,
    H_s,
    H_l
)

Compute the residual ground heat flux that would close the surface energy balance. With all fluxes positive upward (aligned with +z), the energy arriving at the skin from below must balance the radiative and turbulent losses above, so G = R_net + H_s + H_l.

compute_heat_advection_timescale(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Kernel function returning the advective timescale $τ = C Δz / (C_w |j_w|)$ at cell i, j, k, with C the bulk volumetric heat capacity of the cell, C_w that of liquid water, and |j_w| the larger magnitude of the darcy_flux through the two faces of the cell. This is the stability limit of the explicit upwind advection, which changes the cell temperature at the rate C_w j_w / C (m/s). Cells without flow return Inf.

compute_heat_capacity(props, soil)

Compute the bulk heat capacity of the given soil volume.

compute_hydraulic_diffusion_timescale(
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Kernel function returning the hydraulic diffusion timescale $Δz² (∂θ/∂ψ) / K$ at cell i, j, k. The hydraulic conductivity K is evaluated at the cell centre from the local soil composition; the specific moisture capacity ∂θ/∂ψ is the analytic derivative of the soil-water retention curve (from FreezeCurves) evaluated at the matric potential ψₘ, reconstructed from the stored total pressure_head by removing its hydrostatic and elevation components. Cells with zero hydraulic conductivity (fully dry or frozen) impose no restriction and return Inf.

compute_hydraulics!(state, grid, hydrology::SoilHydrology, soil::AbstractSoil, args...)

Compute all state-dependent hydraulic auxiliaries such as hydraulic conductivity and field capacity, and wilting point.

compute_hydraulics!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology,
    args...
)

Kernel function that computes dynamic soil hydraulic properties.

compute_hydraulics!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Computes the unsaturated hydraulic conductivity for RichardsEq configurations of SoilHydrology.

compute_implicit_advection_diagonal(
    i,
    j,
    k,
    grid,
    face_velocity,
    Δt
)

Upwind advection part of the main diagonal of row k, Δt (max(u_{k+½}, 0) - min(u_{k-½}, 0)) / Δz^c_k, i.e. the outflow of ϕ_k through both faces, including the boundary faces of the column. It is nonnegative, which with the nonpositive advective off-diagonals keeps the system an M-matrix.

compute_implicit_advection_lower_diagonal(
    i,
    j,
    k′,
    grid,
    face_velocity,
    Δt
)

Upwind advection part of the lower diagonal at index k′ (row k = k′ + 1), -Δt max(u_{k-½}, 0) / Δz^c_k: an upward velocity through the lower face carries ϕ_{k-1} into cell k. Zero at the bottom (k = 1).

compute_implicit_advection_tendency(
    i,
    j,
    k,
    grid,
    face_velocity,
    ϕ
)

Advective tendency -[∂z(u ϕ_up)]_k, the rate of change of the conserved quantity from the upwind fluxes (compute_implicit_advective_flux) through the faces of cell k, used for the conservative update after each solve. Through the boundary faces of the column only the outflow u ϕ_k is included; the inflow enters through the right-hand side. Zero without advection (face_velocity = nothing).

compute_implicit_advection_upper_diagonal(
    i,
    j,
    k,
    grid,
    face_velocity,
    Δt
)

Upwind advection part of the upper diagonal of row k, Δt min(u_{k+½}, 0) / Δz^c_k: a downward velocity through the upper face carries ϕ_{k+1} into cell k. Zero at the top (k = Nz).

compute_implicit_advective_flux(
    i,
    j,
    k,
    grid,
    face_velocity,
    ϕ
)

First-order upwind advective flux u ϕ_up through face i, j, k, with ϕ_up taken from the cell below the face (k - 1) for u > 0 and above it (k) otherwise, as in the coefficients of the TridiagonalOperator.

compute_implicit_advective_velocity!(face_velocity, i, j, k, grid, fields, energy, hydrology)

Store the advective face velocity C_w j_w (W/m²/K) of the implicit heat operator at face i, j, k in face_velocity, and at the top boundary face for k = Nz, with j_w the water flux (m/s) of the step through the face (compute_advective_water_flux) and C_w the volumetric heat capacity of liquid water, as in the explicit advective_heat_flux. Through the boundary faces, the operator uses only the outflow (see TridiagonalOperator). A no-op without advection (face_velocity = nothing).

compute_implicit_capacity!(
    diagonal_coefficient,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    ::FreeWater,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants
)

Store the regularized apparent heat capacity C_app (J/m³/K) at cell i, j, k in diagonal_coefficient, for the FreeWater freeze curve,

C_app = C + C_pc    if 0 < F < 1 (phase change),
C_app = C           otherwise,

where C is the bulk heat capacity, F the liquid water fraction of the current iterate, and C_pc = 1.5 × 10⁸ J/m³/K a fixed regularization of the latent heat term (a Dirac delta at 0 °C for FreeWater). C_pc only damps the temperature change of phase-change cells between iterates and does not enter the conservative energy update. It is not applied to fully frozen or thawed cells, where it would slow the lagged temperature of cells that keep cooling (warming) and overestimate their conductive fluxes.

compute_implicit_capacity!(
    diagonal_coefficient,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Store the specific moisture capacity C = ∂θ/∂ψ (1/m) of the current iterate at cell i, j, k in diagonal_coefficient (compute_specific_moisture_capacity), bounded below by minimum_moisture_capacity.

compute_implicit_conductivity!(
    face_conductivity,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
)

Store the lagged thermal conductivity at face i, j, k in face_conductivity, interpolated as in the explicit diffusive_heat_flux.

compute_implicit_conductivity!(
    face_conductivity,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})}
) where NF

Store the hydraulic conductivity (m/s) of the current iterate at face i, j, k in face_conductivity, as in the explicit darcy_flux (darcy_conductivity).

compute_implicit_diffusion_diagonal(
    i,
    j,
    k,
    grid,
    face_conductivity,
    face_velocity,
    diagonal_coefficient,
    Δt,
    args...
)

Main diagonal of row k, diag_k = C_k - lower_{k-1} - upper_k + a_k - Δt L_k, with the diffusive off-diagonals lower and upper, the upwind advection part a_k (compute_implicit_advection_diagonal), and L_k from compute_implicit_linear_coefficient evaluated with the trailing args.

compute_implicit_diffusion_lower_diagonal(
    i,
    j,
    k′,
    grid,
    face_conductivity,
    Δt
)

Lower diagonal at index k′, which multiplies ϕ_{k′} in row k = k′ + 1 (the BatchedTridiagonalSolver convention): lower_{k′} = -Δt κ_{k-½} / (Δz^c_k Δz^f_{k-½}), zero at the bottom (k = 1).

compute_implicit_diffusion_tendency(
    i,
    j,
    k,
    grid,
    face_conductivity,
    ϕ
)

Diffusive tendency [∂z(κ ∂z ϕ)]_k, the rate of change of the conserved quantity from the diffusive fluxes through the interior faces of cell k only (consistent with the masked coefficients), used for the conservative update after each solve. Mirrors the explicit diffusive_heat_flux with the lagged face_conductivity.

compute_implicit_diffusion_upper_diagonal(
    i,
    j,
    k,
    grid,
    face_conductivity,
    Δt
)

Upper diagonal of row k, upper_k = -Δt κ_{k+½} / (Δz^c_k Δz^f_{k+½}), zero at the top (k = Nz).

compute_implicit_energy_rhs(
    i,
    j,
    k,
    grid,
    fields,
    diagonal_coefficient,
    explicit_energy
)

Right-hand side C_app Tᵐ - Uᵐ + U⁺ of the linearized system at iterate m of the implicit solver, with C_app = diagonal_coefficient and U⁺ = explicit_energy.

compute_implicit_energy_update(
    i,
    j,
    k,
    grid,
    face_conductivity,
    face_velocity,
    temperature,
    explicit_energy,
    Δt
)

Conservative energy update U = U⁺ + Δt (Q_diff + Q_adv) at i, j, k from the implicit heating rates (W/m³) of conduction, Q_diff = ∂z(κ ∂z T) through the interior faces (compute_implicit_diffusion_tendency), and of advection, Q_adv = -∂z(C_w j_w T_up) through the interior faces and the boundary outflow (compute_implicit_advection_tendency), with U⁺ = explicit_energy, κ = face_conductivity, and C_w j_w = face_velocity (nothing without water flow).

compute_implicit_face_velocity(i, j, k, grid, face_velocity)

Advective face velocity u at face i, j, k of a TridiagonalOperator: face_velocity, or zero(grid) if it is nothing (pure diffusion).

compute_implicit_heat_coefficients!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    fc::FreeWater,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants
)

Store the lagged apparent heat capacity (compute_implicit_capacity!), face conductivity (compute_implicit_conductivity!), and, with water flow, advective face velocity (compute_implicit_advective_velocity!) at i, j, k in out.diagonal_coefficient, out.face_conductivity, and out.face_velocity.

compute_implicit_linear_coefficient(i, j, k, grid, args...)

Coefficient L_k (1/s) of an implicit linear term L ϕ in row k, e.g. a linear sink or a linearized boundary flux; zero by default. args are the trailing arguments of solve! after Δt, for dispatch. A nonzero L must also enter the conservative update as Δt L ϕ.

compute_implicit_meltwater_coefficients!(
    out,
    i,
    j,
    grid,
    fields,
    explicit_water,
    snow::SingleLayerSnow{NF},
    constants::PhysicalConstants,
    Δt
) where NF

Diagonal G′(θᵐ) and right-hand side G′(θᵐ) θᵐ − G(θᵐ) at grid cell i, j of the Newton step for the liquid fraction, with G(θ) = W_i / (1 − θ) + Δt M_r(θ, W_i / (1 − θ)) − W⁺, the storage-limited outflow M_r (compute_meltwater_outflow), the ice content W_i (compute_snow_ice_content), and W⁺ = explicit_water. The linearization point is the liquid fraction 1 − W_i / Wᵐ of the current iterate, bounded to [0, θ⁺]. The coefficients are the identity for a pack without ice or snow, which compute_implicit_snow_water_update! treats separately.

compute_implicit_saturation_update(
    i,
    j,
    k,
    grid,
    fields,
    face_conductivity,
    explicit_saturation,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    Δt
)

Saturation ξ = ξ⁺ + Δt Q_diff / φ at i, j, k, with Q_diff = ∂z(K ∂z Ψ) (1/s) through the interior faces (compute_implicit_diffusion_tendency), ξ⁺ = explicit_saturation, and K = face_conductivity.

compute_implicit_snow_water_update!(
    out,
    i,
    j,
    grid,
    fields,
    explicit_water,
    snow::SingleLayerSnow{NF},
    constants::PhysicalConstants,
    Δt
) where NF

Snow water equivalent W = W_i / (1 − θ) at grid cell i, j from the solved liquid fraction θ = snow_liquid_fraction, bounded to [min(L_c, θ⁺), θ⁺] with the capillary retention L_c, so that W ≤ W⁺ = explicit_water and a Newton step that undershoots the root at the switch between the branches of the outflow drains neither ice nor capillary water, and the drainage meltwater_flux = (W⁺ − W) / Δt. A pack without ice (W_i = 0) takes the exact solution W = max(W⁺ − Δt K_sat, W⁺ / (1 + Δt (1 − L_c) / τ_m)) of its step, and a step without snow (W⁺ ≤ 0) leaves W = W⁺.

compute_implicit_water_flow_coefficients!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Store the capacity (compute_implicit_capacity!) and face conductivity (compute_implicit_conductivity!) at i, j, k in out.

compute_implicit_water_flow_rhs(
    i,
    j,
    k,
    grid,
    fields,
    diagonal_coefficient,
    explicit_saturation,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
)

Right-hand side C Ψᵐ - φ (ξᵐ - ξ⁺) at i, j, k, with C = diagonal_coefficient and ξ⁺ = explicit_saturation.

compute_implicit_water_flux(
    i,
    j,
    k,
    grid,
    water_flux,
    face_conductivity,
    pressure_head
)

Darcy flux q = -K ∂z Ψ (m/s, positive upward) at interior face i, j, k, with K = face_conductivity; at the bottom face, the current water_flux (from tⁿ) is kept.

compute_infiltration(
    runoff::DirectSurfaceRunoff{NF},
    influx,
    sat_top,
    max_infil
) where NF

Compute infiltration from the given influx (water available for infiltration), saturation of the uppermost soil layer sat_top, and the maximum allowed infiltration max_infil.

compute_kinetic_parameters(
    photo::LUEPhotosynthesis{NF},
    T_air::NF
) where NF -> Tuple{Any, Any, Any}

Computes kinetic parameters τ, Kc, Ko based on temperature using Q10 temperature response. Follows enzyme kinetics from [22], Appendix C.

• τ: Rubisco specificity factor (CO₂ to O₂ carboxylation ratio), dimensionless
• Kc: Michaelis-Menten constant for CO₂ [Pa]
• Ko: Michaelis-Menten constant for O₂ [Pa]

References

• [22] Haxeltine & Prentice, Global Biogeochemical Cycles (1996)

compute_latent_heat_flux(
    ::DiagnosedTurbulentFluxes,
    Q_h,
    ρₐ,
    L
)

Compute the latent heat flux as a function of the humidity flux Q_h [m/s], the density ρₐ [kg/m³] of air, and the specific latent heat of vaporization or sublimation L [J/kg].

compute_latent_heat_flux(
    i,
    j,
    grid,
    fields,
    tur::DiagnosedTurbulentFluxes{NF},
    skinT::Terrarium.AbstractSkinTemperature,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
) where NF
compute_latent_heat_flux(
    i,
    j,
    grid,
    fields,
    tur::DiagnosedTurbulentFluxes{NF},
    skinT::Terrarium.AbstractSkinTemperature,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology}
) where NF
compute_latent_heat_flux(
    i,
    j,
    grid,
    fields,
    tur::DiagnosedTurbulentFluxes{NF},
    skinT::Terrarium.AbstractSkinTemperature,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology},
    snow::Union{Nothing, Terrarium.AbstractSnow}
) where NF

Compute the latent heat flux at i, j based on the current skin temperature and atmospheric conditions. When an evapotranspiration scheme is provided, uses the ET-aware humidity flux; otherwise uses bare ground evaporation.

With a snow component, the flux is partitioned by snow-covered area fraction f_snow: the snow-free fraction (1 − f_snow) evaporates from the ground/canopy (latent heat of vaporization), while the snow-covered fraction sublimates from the snowpack (latent heat of sublimation, see compute_snow_sublimation_flux). Without snow it reduces to the bare ground/canopy latent flux.

compute_leaf_area_index(phenol, ϕ, LAI_b)

Compute the instantaneous leaf area index LAI = ϕ·LAI_b from the phenology factor ϕ and the balanced leaf area index LAI_b.

compute_leaf_carbon(vegcarbon_dynamics, traits, LAI_b)

Computes the leaf carbon C_leaf = LAI_b / SLA (kgC/m²) implied by the balanced Leaf Area Index LAI_b and the specific leaf area SLA, [21, Eqs. (76-79)]. Together with compute_root_carbon and compute_stem_carbon, this partitions the single vegetation carbon pool, C_leaf + C_root + C_stem = C_veg; the tissue pools are diagnostic, not prognostic.

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_longwave_up(::DiagnosedRadiativeFluxes, constants::PhysicalConstants, L_down, Ts, ϵ)

Compute outgoing longwave radiation from incoming longwave radiation L_down, surface temperature Ts, and emissivity ϵ.

compute_matric_potential(i, j, k, grid, fields)

Matric potential ψₘ = Ψ - ψ_h - ψ_z (m) at cell i, j, k, reconstructed from the total pressure_head Ψ by removing the hydrostatic head ψ_h = max(0, z₀ - z) below the water table z₀ and the elevation head ψ_z = z - z_ref relative to the soil surface z_ref, as in SoilSaturationPressureClosure.

compute_meltwater_flux!(
    out,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    ::Explicit
)

Store the storage-limited Darcy-type meltwater outflow M_r (m/s SWE) at grid cell i, j in meltwater_flux, evaluated at the diagnosed liquid water fraction and the snow water equivalent of the pack (see compute_meltwater_outflow). A no-op under implicit integration of snow_water_equivalent, where the implicit step sets it.

compute_meltwater_outflow(
    hydraulics::ConstantSnowHydraulics,
    θ_liq::NF,
    W_snow::NF
) where NF

Meltwater outflow M_r (m/s, SWE) from a snowpack with snow water equivalent W_snow (m). Liquid water in excess of the capillary retention L_c drains from the snowpack with a cubic conductivity [5] (in excess-saturation form), K_sat · S*³ with S* = max(θ_liq − L_c, 0) / (1 − L_c), where θ_liq is the liquid fraction of the water substance. This Darcy rate does not depend on the amount of water in the pack, so it is additionally limited by the excess liquid water depth max(θ_liq − L_c, 0) · W_snow drained over the timescale τ_m:

M_r = min(K_sat · S*³, max(θ_liq − L_c, 0) · max(W_snow, 0) / τ_m)

Outflow vanishes smoothly as θ_liq → L_c or W_snow → 0, so a thin pack cannot drain more water than it stores; explicit time steps should not exceed τ_m.

References

• [5] Tarboton, Chowdhury and Jackson (1994)

compute_meltwater_outflow_derivative(
    hydraulics::ConstantSnowHydraulics,
    θ_liq::NF,
    W_ice::NF
) where NF

Derivative dM_r/dθ_liq (m/s) of the meltwater outflow (compute_meltwater_outflow) with respect to the liquid fraction θ_liq at a fixed ice content W_ice (m SWE), i.e. with the snow water equivalent W_snow = W_ice / (1 − θ_liq). It is the derivative of the active branch of the outflow: 3 K_sat S*² / (1 − L_c) where the Darcy rate is the smaller one, and W_ice (1 − L_c) / (τ_m (1 − θ_liq)²) where the storage limit is; zero below the capillary retention L_c. Requires θ_liq < 1.

compute_phenology(i, j, grid, fields, phenol::AbstractPhenology, atmos::AbstractAtmosphere)

Cell-level phenology computation. Implementations return phenology factors and derived LAI at the given index i, j, using atmospheric inputs (e.g. air temperature) where required by the scheme.

compute_phenology!(
    out,
    i,
    j,
    grid,
    fields,
    phenol::PALADYNPhenology,
    atmos::Terrarium.AbstractAtmosphere
)

Mutating wrapper for compute_phenology that stores the result in out.

compute_phenology(
    i,
    j,
    grid,
    fields,
    phenol::PALADYNPhenology,
    atmos::Terrarium.AbstractAtmosphere
) -> Tuple{Any, Any}

Compute the phenology factor and instantaneous leaf area index (LAI) at a single grid point.

compute_phenology_factor(phenol, gdd, T_air)

Compute the phenology factor ϕ [-] from accumulated growing degree days gdd and air temperature T_air.

The cold-deciduous factor combines the growth-phase ramp gdd / gdd_crit [21, Eq. (83)] with a senescence ramp that declines linearly from 1 at T_gdd_base to 0 at T_gdd_base − T_senescence_range. Taking the minimum reproduces all three PALADYN regimes (linear green-up, held at 1 when mature and warm, temperature-driven senescence). The result is then blended with the evergreen value of 1 via the prescribed deciduous fraction.

compute_photosynthesis(i, j, grid, fields, photo::AbstractPhotosynthesis, atmos::AbstractAtmosphere)

Cell-level photosynthesis computation. Implementations compute leaf respiration and net assimilation for a single horizontal cell and return the pair (Rd, An, GPP) or similar outputs as required by the photosynthesis scheme.

compute_photosynthesis!(
    out,
    i,
    j,
    grid,
    fields,
    photo::LUEPhotosynthesis,
    stomcond::Terrarium.AbstractStomatalConductance,
    traits::PlantTraits,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)

Calls compute_photosynthesis and stores the results in out.

compute_photosynthesis(
    i,
    j,
    grid,
    fields,
    photo::LUEPhotosynthesis,
    stomcond::Terrarium.AbstractStomatalConductance,
    traits::PlantTraits,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
) -> Tuple{Any, Any, Any}

Compute photosynthesis, leaf respiration, and gross primary production at a single grid point. Returns instantaneous rates in (kgC/m²/s) for integration by the timestepper.

compute_plant_available_water(
    i, j, k, grid, fields,
    paw::AbstractPlantAvailableWater,
    soil::AbstractSoil
)

Comptue the plant-available water fraction for grid cell i, j and soil layer k.

compute_plant_available_water!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    paw::FieldCapacityLimitedPAW{NF},
    strat::Terrarium.AbstractStratigraphy,
    hydrology::Terrarium.AbstractSoilHydrology,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    args...
) where NF

Mutating wrapper for compute_plant_available_water that stores the result in out.

compute_plant_available_water(
    i,
    j,
    k,
    grid,
    fields,
    paw::FieldCapacityLimitedPAW{NF},
    strat::Terrarium.AbstractStratigraphy,
    hydrology::Terrarium.AbstractSoilHydrology,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Compute the plant avaialble water given the current soil stratigraphy, hydrology, and biogeochemistry state in fields.

compute_precip_ground(
    ::PALADYNCanopyInterception{NF},
    precip,
    I_can,
    R_can
)

Compute rainfall_ground, the rate of rain reaching the ground, following a modified version of [21, Eq. (44)]. Instead of subtracting the tendency, we just directly subtract interception and add the removal rate R_can.

References

• [21] Willeit and Ganopolski, Geoscientific Model Development (2016)

compute_pres_i(photo::LUEPhotosynthesis, λc, pres_a)

Computes intercellular CO2 partial pressure [Pa], [21, Eq. (67)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_quasi_steady_state_temperature(
    i,
    j,
    k,
    grid,
    fields,
    init::QuasiThermalSteadyState
)

Quasi-steady-state temperature T₀ - Qgeo / k_eff * z at cell (i, j, k) from the initial_surface_temperature (°C) and geothermal_heat_flux (W/m²) fields.

compute_radiative_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    rad::PrescribedRadiativeFluxes,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

Compute net radiation and store in auxiliary fields at a grid point.

compute_redistributed_heat!(out, i, j, k, grid, fields, energy, strat, bgc, k_source, Δs)

Change out.internal_energy of cell i, j, k by the sensible heat C_w T φ Δs of a saturation change Δs by redistributed water, with C_w the water heat capacity of the energy process, φ the porosity of cell k, and T the temperature of the cell k_source that the water leaves. A no-op without a SoilThermodynamics energy process.

compute_relative_leaf_area(
    i,
    j,
    grid,
    fields,
    traits::PlantTraits{NF}
) where NF

Return the leaf area index at grid point i, j relative to the maximum leaf area index of the given plant traits, LAI / LAI_max, clamped to the unit interval.

compute_resp10(autoresp)

Computes resp10, the tissue respiration rate at 10 °C per unit tissue nitrogen (kgC/kgN/day).

This is a placeholder that returns the constant parameter resp10. In PALADYN ([21]), the rate of tree PFTs varies with the temperature of the coldest month, between 0.011 for tropical and 0.066 kgC/kgN/day for extratropical trees.

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_respiration!(
    out,
    i,
    j,
    k,
    grid,
    fields,
    soc::SinglePoolSoilCarbon{NF},
    soil::Terrarium.AbstractSoil
) where NF

Compute the respiration rate at index i, j, k and store it in out.respiration_rate.

compute_respiration_assimilation(
    photo::LUEPhotosynthesis{NF},
    traits::PlantTraits{NF},
    constants::MaterialConstants{NF},
    T_air::NF,
    swdown::NF,
    pres::NF,
    co2::NF,
    LAI::NF,
    λc::NF,
    β::NF
) where NF -> Tuple{Any, Any}

Compute and return leaf respiration [kgC/m²/s] and net assimilation [kgC/m²/s] rates.

compute_respiration_rate(
    resp::Terrarium.SoilCarbonRespiration{NF},
    C,
    T
) where NF

Compute the decomposition (heterotrophic respiration) rate (kg/m³/s) of soil organic carbon with density C (kg/m³) at temperature T (°C), R = kref ⋅ Q10^((T - Tref) / 10) ⋅ C.

compute_root_carbon(vegcarbon_dynamics, traits, LAI_b)

Computes the fine root carbon C_root (kgC/m²) implied by the balanced Leaf Area Index LAI_b. Root carbon is assumed equal to leaf carbon (see compute_leaf_carbon), [21, Eqs. (76-79)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_saturation_boundary_flux(
    bc::Oceananigans.BoundaryConditions.BoundaryCondition{<:Flux},
    i,
    j,
    grid,
    clock,
    fields
)

Kernel function returning the flux (positive upward) of the flux boundary condition bc on saturation_water_ice at column i, j, and zero for other boundary conditions.

compute_saturation_closure!(
    state,
    grid,
    closure::SoilSaturationPressureClosure,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    energy::Union{Nothing, Terrarium.AbstractSoilThermodynamics}
) where NF

Saturation–pressure closure of a Richards-equation hydrology: redistribute over- and undersaturation (adjust_saturation_profile!), routing surface excess water to runoff and moving the heat of the redistributed water in the internal_energy of energy (either may be nothing to skip the routing or the heat transfer), then update the water table and compute the pressure_head. Shared by closure! and the implicit Richards step, whose intermediate iterations skip the routing and the heat transfer.

compute_saturation_tendency!(
    saturation_water_ice_tendency,
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration}
)
compute_saturation_tendency!(
    saturation_water_ice_tendency,
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    timestepping::Timestepping
)

compute_sensible_heat_flux(
    i,
    j,
    grid,
    fields,
    tur::DiagnosedTurbulentFluxes,
    skinT::Terrarium.AbstractSkinTemperature,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)

Compute the sensible heat flux at i, j based on the current skin temperature and atmospheric conditions.

compute_sensible_heat_flux(
    ::DiagnosedTurbulentFluxes,
    Q_T,
    ρₐ,
    cₐ
)

Compute the sensible heat flux [W/m²] as a function of the bulk aerodynamic temperature gradient Q_T [K m/s] and the density ρₐ [kg/m³] and specific heat capacity cₐ [J/kg K] of air.

compute_shortwave_up(::DiagnosedRadiativeFluxes, S_down, α)

Compute outgoing shortwave radiation from the incoming shortwave radiation S_down and albedo α.

compute_skin_temperature(
    i,
    j,
    grid,
    fields,
    skinT::ImplicitSkinTemperature{NF},
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    constants::PhysicalConstants,
    snow::Terrarium.AbstractSnow
) where NF

Invert the (linear) area-weighted conduction relation for the implicit skin temperature Ts given the atmosphere-side demanded flux G (see compute_ground_heat_flux_demand), by equating G to the area-weighted sum of the unblended ground and snow-top conductive fluxes, (1 − f_snow)·2κg(Tg − Ts)/Δzg + f_snow·2κsnow(Tsnow − Ts)/dsnow.

compute_skin_temperature(
    i,
    j,
    grid,
    fields,
    skinT::ImplicitSkinTemperature{NF},
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    args...
) where NF

Invert the (linear) ground-only conduction relation for the implicit skin temperature Ts given the atmosphere-side demanded flux G* (see compute_ground_heat_flux_demand): Ts = Tg − G*/(2κg/Δzg). This is the no-snow special case (f_snow = 0) of the snow-aware method below; it is a separate method (rather than a default snow = nothing) purely so it can skip the unused snow_thermal_interface/snow_cover_fraction calls.

compute_skin_temperature_residual!(
    out,
    i,
    j,
    grid,
    fields,
    skinT::ImplicitSkinTemperature,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_skin_temperature_residual!(
    out,
    i,
    j,
    grid,
    fields,
    skinT::ImplicitSkinTemperature,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology}
)
compute_skin_temperature_residual!(
    out,
    i,
    j,
    grid,
    fields,
    skinT::ImplicitSkinTemperature,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology},
    snow::Union{Nothing, Terrarium.AbstractSnow}
)

Surface-energy-balance residual at grid cell i, j, in temperature space: Ts_prev − Ts_implicit, where Ts_implicit is the exact conduction-side inverse (see compute_skin_temperature) of the atmosphere-side demanded flux supplied by the ground heat flux sub-process (see compute_ground_heat_flux_demand), which for DiagnosedGroundHeatFlux is G_demand = R_net(Ts_prev) + H(Ts_prev) + LE(Ts_prev).

compute_snow_basal_heat_flux(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants
)

Snow→soil basal conductive heat flux Q_base [W/m²] at grid cell i, j, positive upward (soil → snow), as the series-resistance conduction between the soil's top half-cell and the snow layer: Q_base = (T_soil − T_snow) / (Δz_soil/(2κ_soil) + d_snow/(2κ_snow)). Both the local soil top-layer thermal conductivity κ_soil and the snow conductivity κ_snow are recovered from their respective composition/density schemes rather than stored. The snow-side conduction thickness is floored at min_snow_conduction_thickness as d_snow → 0, matching the floor used elsewhere in the snow thermodynamics; without the soil-side term this reduces to the previous snow-resistance-only closure, which is only valid once the soil resistance is genuinely negligible next to the (floored) snow's.

compute_snow_cover_fraction(
    cover::FractionalSnowCover{NF},
    swe::NF
) where NF

Sub-grid snow-covered area fraction f_snow = W_snow/(W_snow + W_ref) ∈ [0,1) from the snow water equivalent W_snow [m] and the reference level W_ref (half_coverage), clamping negative SWE (which can occur transiently in the prognostic state) to zero cover.

snow_density(i, j, grid, fields, ::AbstractSnowDensity)

Compute or retrieve the bulk snow density ρ_snow [kg/m³].

compute_snow_depth(
    ::Terrarium.AbstractSnow,
    W_snow::NF,
    ρ_snow::NF,
    ρ_w::NF
) where NF

Snow layer depth d_snow = W_snow·ρ_w/ρ_snow [m], converting the water-equivalent depth W_snow [m] to the physical snow depth using the water density ρ_w and the bulk snow density ρ_snow.

compute_snow_diffusion_timescale(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    constants::PhysicalConstants
)

Kernel function returning the snow thermal diffusion timescale $\max(d_{\text{snow}}, d_{\text{min}})² \, C_{\text{snow}} / κ_{\text{snow}}$ at grid cell i, j, with the bulk thermal conductivity and heat capacity computed from the local snow density and liquid fraction.

compute_snow_energy_tendency(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    atmos::Terrarium.AbstractAtmosphere,
    constants::PhysicalConstants
)

Depth-integrated snow energy tendency [W/m²] at grid cell i, j (all fluxes positive upward):

dŪ_snow/dt = Q_base − Q_top + Q_precip + Q_subl

where Q_top/Q_base are the surface/basal heat fluxes, which like the snow state are grid-cell means scaled by the snow cover fraction f (f·S and f·Q_base, see compute_snow_interface_fluxes!), Q_precip the advected precipitation heat (see compute_snow_precip_heat_flux), and Q_subl an advective correction for sublimation.

The sublimation correction Q_subl = ρ_w·L_sl·E_subl is required because the latent heat flux carries the full sublimation enthalpy ρ_w·L_sg·E_subl, whereas the mass leaving the snowpack departs as ice, whose specific enthalpy relative to the liquid-water reference is −L_sl. Adding back ρ_w·L_sl·E_subl leaves the snowpack with a net loss of ρ_w·(L_sg − L_sl)·E_subl = ρ_w·L_lg·E_subl, the vaporization enthalpy carried by the departing vapor.

Note that no explicit meltwater energy term appears because meltwater drains as liquid water at 0 °C, which is the zero-enthalpy reference (U = 0) of the FreeWater closure, so it carries no enthalpy out of the snowpack.

The conductive/sublimation terms Q_base - Q_top + Q_subl are gated to zero when the pack is empty at the start of the step (W == 0), since they are otherwise evaluated against a degenerate zero-depth layer (Q_top/Q_base conduction targets floored at min_conduction_thickness, giving an arbitrarily large flux relative to a genuinely-thin pack's tiny heat capacity). Q_prcp is not gated: it is a pure source term from the precipitation rate alone, well-defined regardless of the existing pack depth. Gating it too — as a previous implementation did, multiplying the whole dUdt including Q_prcp by W > 0 — zeroed the energy input on the very first step any snow accumulates onto bare ground (W is read before that step's mass increment), leaving the freshly-fallen, physically cold snow with zero recorded enthalpy. Under the FreeWater closure (U = 0 ≡ liquid water at 0 °C) that reads as instantaneously fully melted, triggering spurious immediate Darcy drainage regardless of how cold the air actually was.

compute_snow_ice_content(
    i,
    j,
    grid,
    fields,
    explicit_water,
    constants::PhysicalConstants
)

Ice content W_i = −Ū⁺ / (ρ_w L_f) (m SWE) of the snowpack at grid cell i, j, clamped to [0, W⁺], from the explicit stage W⁺ = explicit_water and the snow energy Ū⁺; unchanged by the meltwater drainage.

compute_snow_interface_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)

Diagnose the snow↔surface/soil coupling fluxes at grid cell i, j, run after the surface energy balance: the blended soil-top heat flux f·Q_base + (1 − f)·G (with G the bare-ground ground_heat_flux; used only as the soil's top boundary condition), the grid-cell-mean snow-base flux f·Q_base (see compute_snow_basal_heat_flux), the grid-cell-mean snow-top conductive flux f·S (see compute_snow_surface_heat_flux), and the snow surface sublimation rate (see compute_snow_sublimation_flux, the same snow-fraction vapor flux the surface energy balance uses for the latent-flux partition). All snowpack fluxes are grid-cell means scaled by the snow cover fraction f (as is the sublimation) since the snow state is a grid-cell mean. The column (soil + snow) then loses (1 − f)·G + f·S, the flux delivered to the skin.

compute_snow_interface_fluxes!(
    state,
    grid,
    ::Nothing,
    args...
)

Launch compute_snow_interface_fluxes! to diagnose the snow↔surface/soil coupling fluxes from the snow state and the surface energy balance outputs: the blended soil-top heat flux (soil_heat_flux), the grid-cell-mean snow-base flux (snow_basal_heat_flux), the grid-cell-mean snow-top conductive flux (snow_surface_heat_flux; both drive the snowpack's own energy tendency), and the snow surface sublimation rate (sublimation, the snow-fraction bulk-aerodynamic vapor flux at the converged skin temperature). Must run after the surface energy balance (which sets ground_heat_flux and the skin temperature). No-op when there is no snowpack (snow === nothing).

compute_snow_precip_heat_flux(
    ::Terrarium.AbstractSnow,
    constants::PhysicalConstants,
    P_s::NF,
    R_on_snow::NF,
    T_air::NF
) where NF

Advected heat flux [W/m²] carried into the snowpack by precipitation, relative to liquid water at 0 °C (the U = 0 reference of the FreeWater enthalpy closure). Fresh snow P_s arrives as ice, which sits L_sl below the liquid reference, plus sensible heat for T_air < 0; rain-on-snow R_on_snow arrives as liquid carrying only its sensible heat for T_air > 0. The latent heat released when rain refreezes in a cold pack is captured implicitly by the enthalpy closure, so it is not added here (adding L_sl would double-count relative to the liquid-water reference).

compute_snow_properties!(
    out,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    constants::PhysicalConstants
)
compute_snow_properties!(
    out,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    water_timestepping::Timestepping
)

Compute the snow depth and cover fraction at grid cell i, j, and the meltwater outflow if snow_water_equivalent is integrated explicitly (see compute_meltwater_flux!).

compute_snow_sublimation_flux(
    i,
    j,
    grid,
    fields,
    ::Nothing,
    atmos,
    constants,
    skinT
)

Snow-surface sublimation rate [m/s SWE] at grid cell i, j. The snow surface is treated as saturated: a bulk-aerodynamic vapor flux Δq/rₐ evaluated at the skin temperature, with Δq taken over ice for a sub-freezing surface (the saturation humidity already dispatches over ice for T ≤ 0 — see saturation_specific_humidity_vapor). The water-vapor mass flux ρₐ·Δq/rₐ is converted to a snow-water-equivalent rate via ρ_w. Zero without snow (i.e. when snow === nothing).

compute_snow_surface_heat_flux(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    constants::PhysicalConstants
)

Snow-top conductive heat flux [W/m²] at grid cell i, j (positive upward, i.e. energy leaving the snow top into the skin), S = 2κ_snow(T_snow − T_skin)/max(d_snow, d_min) — the same unblended snow conduction target the skin-temperature solve closes against (see Terrarium.snow_thermal_interface). This is a per-unit-snow-area flux; the snowpack energy tendency (compute_snow_energy_tendency's Q_top) uses the grid-cell mean f_snow·S, which is applied at the call site in compute_snow_interface_fluxes!.

compute_snow_tendencies!(
    tendencies,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    atmos::Terrarium.AbstractAtmosphere,
    constants::PhysicalConstants
)
compute_snow_tendencies!(
    tendencies,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    atmos::Terrarium.AbstractAtmosphere,
    constants::PhysicalConstants,
    water_timestepping::Timestepping
)

Accumulate the snow water-equivalent and depth-integrated energy tendencies at grid cell i, j from the mass and energy balances (see compute_snow_water_tendency and compute_snow_energy_tendency).

compute_snow_volumetric_energy(
    Ū_snow::NF,
    d_snow::NF,
    d_min::NF
) where NF

Volumetric snow internal energy U_snow = Ū_snow/max(d_snow, d_min) [J/m³] from the depth-integrated energy Ū_snow [J/m²] and the snow depth d_snow [m]. The depth is floored at the minimum thermal thickness d_min (see min_snow_conduction_thickness) rather than a machine-eps offset. This bounds the snow temperature recovered downstream: with only an eps offset, any residual Ū_snow over a vanishing d_snow gives a huge U_snow and hence a snow temperature far below physical bounds, which then corrupts the (cover-fraction-blended) skin temperature and basal heat flux. Flooring at d_min gives a thin snowpack a bounded effective heat capacity C_snow·d_min, so it stores negligible energy and is thermally transient — the ground heat flux passes essentially unmediated to the soil (f_snow → 0 in the blend) while the pack still stores mass and modifies the surface albedo.

compute_snow_volumetric_heat_capacity(
    ::Terrarium.AbstractSnow{NF},
    constants::PhysicalConstants,
    ρ_snow::NF,
    liq::NF
) where NF

Bulk volumetric heat capacity C_snow [J/m³/K] of the snowpack treated as an ice–liquid–air mixture, given the bulk snow density ρ_snow [kg/m³] and the liquid water fraction liq ∈ [0,1] of the water substance. The water-substance mass per unit snow volume is ρ_snow, of which (1 − liq) is ice and liq is liquid; the remaining void space is dry air. The corresponding constituent volume fractions are

\[\begin{aligned} θ_{ice} &= ρ_{snow}·(1 − liq)/ρ_i,\ θ_{liq} &= ρ_{snow}·liq/ρ_w,\ θ_{air} &= 1 − θ_{ice} − θ_{liq} \end{aligned}\]

and the heat capacity is the volume-weighted sum over the ice and liquid constituents, C_snow = cp_i·ρ_i·θ_ice + cp_w·ρ_w·θ_liq. The air's own sensible-heat storage (≈0.1%) is neglected; note the air is still reflected in the bulk density ρ_snow < ρ_ice.

compute_snow_water_tendency(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow{NF},
    atmos::Terrarium.AbstractAtmosphere
) where NF
compute_snow_water_tendency(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow{NF},
    atmos::Terrarium.AbstractAtmosphere,
    water_timestepping::Timestepping
) where NF

Snow water equivalent (SWE) tendency (m/s) at grid cell i, j:

dW_snow/dt = S + R_snow − M − E_subl

where S is snowfall, R_snow = f_snow · rainfall the rain intercepted by the snow-covered fraction, M the Darcy meltwater outflow (see snow_meltwater_flux), and E_subl the sublimation rate. M is omitted under implicit integration of snow_water_equivalent, where the implicit step drains the meltwater.

compute_soc_advective_flux(
    i,
    j,
    k,
    grid,
    fields,
    transport::Terrarium.SoilCarbonTransport
)

Compute the advective soil organic carbon flux q_a = ω ⋅ C (kg/m²/s) at index i, j, k.

compute_soc_diffusive_flux(
    i,
    j,
    k,
    grid,
    fields,
    transport::Terrarium.SoilCarbonTransport
)

Compute the diffusive soil organic carbon flux qd = -Db ⋅ ∂C/∂z (kg/m²/s) at the lower cell face of index i, j, k.

compute_soc_tendency!(
    tend,
    i,
    j,
    k,
    grid,
    fields,
    soc::SinglePoolSoilCarbon,
    soil::Terrarium.AbstractSoil
)

Compute the soil organic carbon tendency at index i, j, k and add it to tend.density_soc.

compute_soc_tendency(
    i,
    j,
    k,
    grid,
    fields,
    soc::SinglePoolSoilCarbon,
    soil::Terrarium.AbstractSoil
)

Compute the soil organic carbon tendency (kg/m³/s) at index i, j, k, ∂C/∂t = -∂qd/∂z - ∂qa/∂z - R, where qd and qa are the diffusive and advective fluxes and R is the respiration rate. Litter inputs enter through the top boundary condition.

compute_specific_humidity_difference(
    i,
    j,
    grid,
    fields,
    atmos::Terrarium.AbstractAtmosphere,
    c::PhysicalConstants,
    T
)

Computes the specific humidity difference [kg/kg] between a saturated surface at temperature T [°C] and the current atmospheric fields.

compute_specific_moisture_capacity(
    i,
    j,
    k,
    grid,
    fields,
    hydrology,
    strat,
    bgc
)

Specific moisture capacity ∂θ/∂ψ (1/m) at cell i, j, k: the derivative of the soil water retention curve (specific_moisture_capacity) at the matric potential (compute_matric_potential) reconstructed from the total pressure_head. It vanishes in saturated cells.

compute_stem_carbon(vegcarbon_dynamics, traits, LAI_b)

Computes the stem carbon C_stem = awl · LAI_b (kgC/m²) implied by the balanced Leaf Area Index LAI_b and the allometric coefficient awl (assuming bwl = 1), [21, Eqs. (76-79)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_stomatal_conductance(
    i, j, grid, fields,
    stomcond::AbstractStomatalConductance,
    photo::AbstractPhotosynthesis,
    atmos::AbstractAtmosphere,
    constants::PhysicalConstants,
    args...
)

Cell-level stomatal conductance computation. Returns stomatal/canopy conductance and internal CO₂ ratio for the specified cell.

compute_stomatal_conductance!(
    out,
    i,
    j,
    grid,
    fields,
    stomcond::MedlynStomatalConductance{NF},
    traits::PlantTraits{NF},
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
) where NF

Calls compute_stomatal_conductance and stores the result in out.

compute_stomatal_conductance(
    i,
    j,
    grid,
    fields,
    stomcond::MedlynStomatalConductance{NF},
    traits::PlantTraits{NF},
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
) where NF

Compute stomatal conductance (gstm) and leaf-to-air CO₂ ratio (λc) at a grid point. Returns tuple (gstm, λc) for use in photosynthesis and transpiration calculations.

compute_stomatal_conductance(
    stomcond::MedlynStomatalConductance{NF},
    traits::PlantTraits{NF},
    constants::PhysicalConstants{NF},
    vpd,
    T_air,
    pres,
    co2_ppm,
    An,
    LAI,
    β
) where NF

Compute canopy-level water conductance (m/s) from the [24] optimal stomatal conductance model, gₛ = g₀ + D (1 + g₁/√VPD) Aₙ/cₐ, where the diffusivity ratio D is that of the corrigendum [25]. Includes minimum conductance and light extinction effects based on LAI, scaled by soil moisture factor β.

References

• [24] Medlyn et al., Global Change Biology (2011)
• [25] Medlyn et al., Global Change Biology (2012)

compute_surface_drainage(
    runoff::DirectSurfaceRunoff{NF},
    surface_excess_water
) where NF

Compute surface drainage flux from the current surface_excess_water resevoir state.

compute_surface_energy_fluxes!(
    state,
    grid,
    seb::SurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_surface_energy_fluxes!(
    state,
    grid,
    seb::SurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology},
    args...
)

Compute the surface energy fluxes on grid based on the current atmospheric state.

compute_surface_energy_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    seb::SurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_surface_energy_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    seb::SurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology}
)
compute_surface_energy_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    seb::SurfaceEnergyBalance,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology},
    snow::Union{Nothing, Terrarium.AbstractSnow},
    args...
)

Fused kernel function that computes the radiative and turbulent fluxes, as well as the ground heat flux based on the current skin temperature and humidity fluxes.

compute_surface_energy_fluxes!(state, grid, ::AbstractSurfaceEnergyBalance, args...)

Compute the surface energy fluxes and skin temperature from the current state and grid. The required args are implementation dependent.

compute_surface_excess_water_tendency!(
    tendencies,
    i,
    j,
    grid,
    fields,
    runoff::DirectSurfaceRunoff{NF}
) where NF

Kernel function for computing the tendency of the prognostic surface_excess_water pool in the grid column i, j. The pool is drained at the surface drainage rate D = compute_surface_drainage, capped so that a single step cannot remove more water than is present (min(D, S)). The tendency is returned as a negative rate (a removal) consistent with the positive-upward flux convention used by the other surface hydrology tendencies, so that surface_excess_water += ∂S∂t * Δt draws the pool down.

compute_surface_humidity_flux(
    i,
    j,
    grid,
    fields,
    evaporation::BareGroundEvaporation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

Evaluate the bare-ground, kinematic surface humidity flux [m/s] from the current skin temperature and evaporation conductance in fields.

compute_surface_humidity_flux(
    i,
    j,
    grid,
    fields,
    evtr::PALADYNCanopyEvapotranspiration,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

Compute the kinematic surface humidity flux [m/s] from the current skin temperature and conductances in fields.

Compute unscaled transpiration, evaporation_ground, and evaporation_canopy fluxes on grid. Following the implementation of PALADYN, Qh_can is clamped to be strictly positive (no canopy dew formation).

Compute the net radiation budget given incoming and outgoing shortwave and longwave radiation.

compute_surface_runoff!(
    out, i, j, grid, fields,
    runoff::DirectSurfaceRunoff{NF},
    canopy_interception::AbstractCanopyInterception,
    soil_hydrology::AbstractSoilHydrology
) where {NF}

Compute surface runoff in grid cell i, j and store the result in out.surface_runoff. Runoff is computed based on the ground-reaching precipitation rate provided by canopy_interception and the soil saturation state provided by soil_hydrology.

compute_surface_runoff(
    runoff::DirectSurfaceRunoff,
    influx,
    surface_drainage,
    infil
)

Compute surface runoff as precipitation + surface_drainage - infiltration.

compute_surface_upwelling_radiation(
    i,
    j,
    grid,
    fields,
    rad::PrescribedRadiativeFluxes,
    args...
) -> NamedTuple{(:surface_shortwave_up, :surface_longwave_up), <:Tuple{Any, Any}}

Compute upwelling shortwave and longwave radiation at a grid point.

compute_temperature_stress(
    photo::LUEPhotosynthesis{NF},
    T_air::NF
) where NF

Computes the temperature stress factor T_stress based on the air temperature.

compute_tendencies!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants
)
compute_tendencies!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology}
)
compute_tendencies!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology},
    timestepper::Union{Nothing, Terrarium.AbstractTimeStepper}
)

Compute tendencies for soil energy, water, and carbon state variables on grid based on the given values in constants.

An optional surface_hydrology process may be supplied (by the coupled LandModel) so that its evapotranspiration scheme (get_evapotranspiration) is applied as a sink term in the soil water tendency. Without it (standalone soil), no evapotranspiration is removed from the soil.

The parent model's timestepper resolves the Timestepping classes of internal_energy and saturation_water_ice, which are passed on to the energy and hydrology compute_tendencies!. Without a timestepper (nothing) both are Explicit().

compute_tendencies!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil
)
compute_tendencies!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::Union{Nothing, PhysicalConstants}
)
compute_tendencies!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::Union{Nothing, PhysicalConstants},
    timestepping::Timestepping
)

Compute the explicitly integrated part of the internal energy tendency of energy. The Timestepping argument is the class of internal_energy (see timestepping(state, ::Val)), i.e. how the variable is integrated, not which part of the tendency is computed. Explicit: the whole tendency is integrated explicitly, i.e. the full divergence of the diffusive and advective heat fluxes. Implicit: the fluxes through the interior faces are integrated by the implicit sub-stepper (BackwardEuler), so only the diffusive fluxes and the advective inflow through the boundary faces remain. Flux boundary conditions are included in both cases.

compute_tendencies!(state, model::AbstractModel)

Compute tendencies for all prognostic state variables for model stored in the given state. This method should be called after compute_auxiliary!.

compute_tendencies!(state, grid, process::AbstractProcess, args...)

Compute the tendencies of all prognostic state variables for the given process on grid. Implementations of AbstractProcess may define additional args that correspond to different process coupling interfaces.

compute_tendencies!(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)
compute_tendencies!(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    timestepper::Union{Nothing, Terrarium.AbstractTimeStepper}
)

Compute the snow water-equivalent and depth-integrated energy tendencies (see compute_snow_tendencies!), driven by the atmospheric precipitation inputs and the prescribed surface/basal heat fluxes and sublimation. The meltwater outflow is omitted if snow_water_equivalent is integrated implicitly by timestepper.

compute_tendencies!(
    state,
    grid,
    runoff::DirectSurfaceRunoff,
    args...
)

compute_tendencies!(
    state,
    grid,
    canopy_interception::PALADYNCanopyInterception,
    evapotranspiration::Terrarium.AbstractEvapotranspiration,
    args...
)

compute_tendencies!(
    state,
    grid,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits,
    args...
)

compute_tendencies!(
    state,
    grid,
    phenol::PALADYNPhenology,
    atmos::Terrarium.AbstractAtmosphere
)

compute_tendencies!(
    state,
    grid,
    veg_dynamics::PALADYNVegetationDynamics,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits,
    args...
)

compute_tendencies!(
    state,
    grid,
    veg::PrognosticVegetation,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere,
    args...
)

Compute tendencies for carbon and vegetation dynamics.

compute_tendencies!(
    state,
    grid,
    soc::SinglePoolSoilCarbon,
    soil::Terrarium.AbstractSoil,
    args...
)

Compute the tendency of the soil organic carbon density in each soil layer.

compute_tendencies!(
    state,
    grid,
    hydrology::SoilHydrology,
    soil::Terrarium.AbstractSoil,
    args...
)

compute_tendencies!(
    state,
    grid,
    hydrology::SurfaceHydrology,
    args...
)

compute_tendencies!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants
) where NF
compute_tendencies!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration}
) where NF
compute_tendencies!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    timestepping::Timestepping
) where NF

Compute the saturation tendency of the Richards-equation hydrology, with the evapotranspiration sink of evtr (if given) and the user-defined vwc_forcing as source terms. timestepping is the resolved Timestepping class of saturation_water_ice (see timestepping(state, ::Val)). Under Explicit the full divergence of the Darcy flux is added. Under Implicit the fluxes through the interior faces are omitted (integrated by the implicit sub-stepper, see BackwardEuler), leaving the fluxes through the two boundary faces (see compute_boundary_darcy_flux_divergence). The flux boundary conditions on saturation_water_ice (e.g. InfiltrationFlux) from compute_boundary_conditions! and the source terms are part of the tendency in both cases.

The Darcy fluxes at tⁿ through all faces are stored in water_flux (compute_water_flux), with which the soil energy advects heat. Under Implicit time stepping, the implicit step overwrites the interior fluxes with those of its conservative update (see implicit_update!).

compute_thermal_conductivity(i, j, k, grid, ::AbstractThermodynamics, args...)

Compute the thermal conductivity at index i, j, k.

compute_thermal_conductivity(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
)

compute_thermal_conductivity(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

compute_thermal_conductivity(
    snow::SingleLayerSnow,
    constants::MaterialConstants,
    ρ_snow
)

Bulk snow thermal conductivity κ_snow [W/m/K], delegating to the process's thermal conductivity scheme with the bulk density ρ_snow.

compute_thermal_conductivity(props, soil)

Compute the bulk thermal conductivity of the given soil volume.

compute_thermal_conductivity(
    cond::PowerLawSnowThermalConductivity{NF},
    constants::MaterialConstants{NF},
    ρ_snow::NF
) where NF

Bulk snow thermal conductivity via the density power law κ_snow = a·(ρ_snow/ρ_w)^b ([18]).

compute_thermal_diffusion_timescale(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
)

Kernel function returning the thermal diffusion timescale $Δz² C / κ$ at cell i, j, k, with the bulk thermal conductivity κ and heat capacity C computed from the local soil composition.

compute_turbulent_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    tur::DiagnosedTurbulentFluxes,
    skinT,
    constants,
    atmos,
    hydrology,
    snow
)

Compute the turbulent (sensible and latent) heat fluxes from the current skin temperature and store them into the auxiliary output fields out.

compute_turbulent_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    ::PrescribedTurbulentFluxes,
    args...
)

Prescribed turbulent fluxes are input fields supplied by an external coupler (already present in fields), so there is nothing to compute or store.

compute_vadose_zone_saturation(
    i,
    j,
    k,
    grid,
    fields,
    saturation
)

Saturation in the vadose zone at cell (i, j, k) from the vadose_zone_saturation input field. For a Saturation or a plain value, the field holds the saturation and is returned directly; for a VolumetricWaterContent, it holds the volumetric water content θ, which is divided by the soil porosity ϕ and capped at full saturation.

compute_vapor_pressure_deficit(
    i,
    j,
    grid,
    fields,
    atmos::Terrarium.AbstractAtmosphere,
    c::PhysicalConstants
)

Computes the vapor pressure deficit (VPD) [Pa] at atmospheric reference level given the current atmospheric fields

compute_vapor_pressure_difference(
    i,
    j,
    grid,
    fields,
    atmos::Terrarium.AbstractAtmosphere,
    c::PhysicalConstants,
    T
)

Computes the vapor pressure difference [Pa] between a saturated surface at temperature T [°C] and the current atmospheric fields.

compute_veg_carbon_auxiliary!(
    out,
    i,
    j,
    grid,
    fields,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits
)

Mutating wrapper for compute_balanced_leaf_area_index that stores the result in out.balanced_leaf_area_index.

compute_veg_carbon_tendencies!(
    tend,
    i,
    j,
    grid,
    fields,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits
)

Calls compute_veg_carbon_tendency and stores the result in out.

compute_veg_carbon_tendency(i, j, grid, fields, vegcarbon::AbstractVegetationCarbonDynamics)

Cell-level vegetation-carbon tendency computation. Implementations compute the tendency for the total vegetation carbon pool at the given index i, j.

compute_veg_carbon_tendency(
    i,
    j,
    grid,
    fields,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits
)

Compute the tendency for the carbon vegetation pool given fields LAI_b and NPP.

compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology,
    constants,
    evapotranspiration
)
compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology,
    constants,
    evapotranspiration,
    timestepping
)

Compute the volumetric water content (VWC) tendency at grid cell i, j k according to the Richardson-Richards equation. Note that the VWC tendency is not scaled by the porosity and is thus not the same as the saturation tendency. Like the energy tendency, it is the explicitly integrated part of the tendency, and timestepping is the class of saturation_water_ice: under Explicit (the default) the full divergence of the Darcy flux is included; under Implicit only the boundary-face fluxes are (compute_boundary_darcy_flux_divergence), the interior ones being integrated by the implicit sub-stepper.

compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration}
) where NF
compute_volumetric_water_content_tendency(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    ::Timestepping
) where NF

Compute the volumetric water content (VWC) tendency at grid cell i, j k f. Note that the VWC tendency is not scaled by the porosity and is thus not the same as the saturation tendency.

compute_water_flux(i, j, k, grid, ψ, K)

Darcy flux q = -K ∂Ψ/∂z (m/s, positive upward) through face i, j, k (darcy_flux), stored in water_flux by compute_tendencies! at tⁿ. At the boundary faces it is the flux implied by the boundary conditions on pressure_head (zero for no-flux, see compute_boundary_darcy_flux_divergence).

compute_water_outflow_timescale(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry,
    constants::PhysicalConstants,
    evtr::Union{Nothing, Terrarium.AbstractEvapotranspiration},
    timestepping::Timestepping,
    saturation_bcs::NamedTuple
) where NF

Kernel function returning the explicit water outflow timescale $τ = φ (s - s_res) Δz / q_out$ at cell i, j, k, with porosity φ, saturation s, the residual_saturation s_res, and the outflow $q_out = -Δz G_θ$ (m/s) where the explicit volumetric water content tendency $G_θ$ (compute_volumetric_water_content_tendency plus the flux boundary conditions saturation_bcs on saturation_water_ice, see compute_saturation_boundary_flux) is negative. Cells without net explicit outflow return Inf.

compute_water_table!(
    water_table,
    i,
    j,
    grid,
    sat,
    ::SoilHydrology{NF, VerticalFlow, SaturationClosure, SoilHydraulics} where {VerticalFlow<:Terrarium.AbstractVerticalFlow, SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})}
) where NF

Kernel function that diagnoses the water table at grid cell i, j given the current soil saturation profile.

compute_water_table_saturation(
    i,
    j,
    k,
    grid,
    fields,
    init::SaturationWaterTable,
    args...
)

Initial saturation at cell (i, j, k): unity at and below the water_table_depth, otherwise the vadose zone saturation (see compute_vadose_zone_saturation).

compute_Γ_star(
    photo::LUEPhotosynthesis{NF},
    τ::NF,
    pres_O2::NF
) where NF

Computes the CO₂ compensation point Γ_star [Pa]. The intercellular CO₂ partial pressure at which gross photosynthesis equals respiration. Follows [21, Eq. (C6)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_Λ_loc(vegcarbon_dynamics, traits, LAI_b)

Computes the local litterfall rate Λ_loc based on the balanced Leaf Area Index LAI_b (assuming evergreen PFTs), [21, Eq. (75)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_γv(veg_dynamics)

Computes the disturbance rateγv, [21, Eq. (80)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_λ_NPP(vegcarbon_dynamics, traits, LAI_b)

Computes λ_NPP,a factor determining the partitioning of NPP between increase of vegetation carbon of the existing vegetated area and spreading of the given PFT based on the balanced Leaf Area Index LAI_b, [21, Eq. (74)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

compute_λc(stomcond, vpd)

Computes the ratio of leaf-internal and air CO2 concentration λc = cᵢ/cₐ, derived from the optimal stomatal conductance model ([24]) in the corrected form of [25].

References

• [24] Medlyn et al., Global Change Biology (2011)
• [25] Medlyn et al., Global Change Biology (2012)

compute_ν_star(veg_dynamics, ν)

Computes ν_star which is the maximum of the current vegetation fraction ν and the seed fraction ν_seed [-], to ensure that a PFT is always seeded.

compute_ν_tendencies!(
    tend,
    i,
    j,
    grid,
    fields,
    veg_dynamics::PALADYNVegetationDynamics,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits
)

Mutating wrapper for compute_ν_tendency that stores the result in tend.

compute_ν_tendency(
    i, j, grid, fields,
    veg_dynamics::AbstractVegetationDynamics,
    vegcarbon::AbstractVegetationCarbonDynamics
)

Cell-level vegetation-fraction tendency computation used by vegetation dynamics. Implementations compute the local tendency value for ν at the given index i, j.

compute_ν_tendency(
    i,
    j,
    grid,
    fields,
    veg_dynamics::PALADYNVegetationDynamics,
    vegcarbon_dynamics::PALADYNCarbonDynamics,
    traits::PlantTraits
)

Compute vegetation area fraction tendency at a single grid point from NPP-productivity and disturbance rates.

compute_ν_tendency(
    veg_dynamics,
    vegcarbon_dynamics,
    traits,
    LAI_b,
    C_veg,
    NPP,
    ν
)

Computes the vegetation fraction tendency for a single PFT, [21, Eq. (73)].

References

• [21] Willeit & Ganopolski, Geoscientific Model Development (2016)

convert_dt(Δt)

Convert Δts to the given target Period or numeric type. Numeric types are always assumed to be in seconds.

coordinate_eltype(z)

Return the number format to use for a grid with vertical coordinate z when none is given explicitly. Oceananigans discretizations compute their cell interfaces in Float64, so only a coordinate given directly as an array of interfaces carries a number format of its own.

create_land_grid(grid::AbstractGrid, soil, snow, vegetation)::AbstractLandGrid

Construct the LandGrid for a model discretized on the spatial discretization grid whose ground, snow, and canopy domains are determined by the soil, snow, and vegetation process components. This is the standard interface through which model constructors turn the grid given by the user into the land grid stored by the model.

The default implementation returns a ground-only land grid, i.e. one for which the snow and canopy domains are not vertically resolved. This is correct for all currently implemented snow and vegetation schemes, which are 0D; schemes which resolve a vertical snow or canopy profile should add a method here which builds the corresponding domain grids.

darcy_conductivity(i, j, k, grid, ψ, K)

Hydraulic conductivity (m/s) used by the darcy_flux at face k: the minimum of K (defined on faces) at face k and at the adjacent face in the direction of the gradient of the total pressure head ψ (face k - 1 for ∂ψ/∂z < 0, else k + 1).

darcy_flux(i, j, k, grid, ψ, K)

Kernel function for computing the Darcy flux q = -K ∂ψ/∂z (m/s, positive upward) over layer faces from the total pressure head ψ and the face hydraulic conductivity K (see darcy_conductivity).

daytime_length(i, j, grid, fields, ::AbstractAtmosphere{PR, <:LongShortWaveRadiation})

Retrieve the length of the day (in hours) at grid cell i, j. Defaults to a constant 12 hours if no input is provided.

debughook!(args...)

Provides a "hook" for handling debug calls from relevant callsites. Default implementations for Field and NamedTuple (assumed to be of Fields) simply forward to checkfinite!.

debugsite!(args...)

Utility method that forwards args to debughook! if and only if debug mode is enabled. Debug mode is set by the global variable DEBUG which can be toggled by the user facing API debug!.

deduplicate_vars(vars::Tuple{Vararg{Union{AbstractVariable, Namespace}}})

Type-stable equivalent of deduplicate for tuples of AbstractVariables and Namespaces.

default_clock(
    model::Terrarium.AbstractModel{NF}
) where NF -> Clock{<:Any, <:Any, Float64, Int64, Int64}

Return the default Clock used by initialize for the given model. The generic method returns a plain host clock starting at time zero; architecture extensions may specialize on the model's grid to return an architecture-specific clock (e.g. TerrariumReactantExt returns a traced ConcreteRNumber-backed clock so that time advances inside the compiled step).

default_dt(timestepper::AbstractTimeStepper)

Get the current timestep size for the time stepper.

default_dt(integrator::ModelIntegrator)

Return the default timestep size for the given integrator, taken from its model's timestepper.

default_skin_temperature_solver(
    ::Type{NF}
) where NF -> Terrarium.RootSolver{NF, RootSolvers.NewtonsMethod{T}, RootSolvers.CompactSolution, RootSolvers.ResidualTolerance{NF}} where {NF, T<:Union{Real, AbstractArray}}

Construct the default solver for the implicit skin temperature: a Newton root-finder (RootSolver backed by RootSolvers.jl) with a small iteration budget.

default_timestepper(::Type{NF})
default_timestepper(::Type{<:AbstractModel}, ::Type{NF})

Default timestepper for models: a single explicit ForwardEuler time stepper. Models may specialize it on their type: LandModel defaults to IMEX(ForwardEuler(NF; Δt = 900), BackwardEuler(NF)), and SoilModel to IMEX(ForwardEuler(NF; Δt = 3600), BackwardEuler(NF)) if its soil water flow follows the Richards equation and to ForwardEuler(NF) otherwise.

default_timestepping(hydrology::AbstractSoilHydrology)

Default Timestepping class of saturation_water_ice under an AbstractIMEX timestepper: Implicit for the Richards equation (SoilHydrology{NF, RichardsEq}) and Explicit otherwise.

julia> Terrarium.default_timestepping(SoilHydrology(Float64, RichardsEq()))
Implicit()

julia> Terrarium.default_timestepping(SoilHydrology(Float64, NoFlow()))
Explicit()

default_unsat_hydraulic_cond(
    ::Type{NF},
    ::SWRC
) where NF -> UnsatKLinear

Return the unsaturated hydraulic conductivity formulation that matches the soil water retention curve swrc: UnsatKBrooksCorey for BrooksCorey, UnsatKVanGenuchten for VanGenuchten, and UnsatKLinear otherwise.

julia> Terrarium.default_unsat_hydraulic_cond(Float32, BrooksCorey())
UnsatKBrooksCorey{Float32}(7.0f0)

julia> Terrarium.default_unsat_hydraulic_cond(Float32, VanGenuchten())
UnsatKVanGenuchten{Float32}(7.0f0)

density_pure_soc(bgc::AbstractSoilBiogeochemistry)

Return the assumed constant density of pure organic material. The default implementation assumes there to be a property ρ_org defined on the type of bgc.

density_soc(i, j, k, grid, fields, bgc::AbstractSoilBiogeochemistry{NF}) where {NF}

Return the bulk soil organic carbon density ρsoc (kg/m³) at index i, j, k, i.e. the mass of organic material per unit volume of soil. The default implementation returns zero, corresponding to a purely mineral soil. Implementations of AbstractSoilBiogeochemistry should override this method; see `organicfraction` for how ρ_soc is converted to the organic fraction of the soil matrix.

describe_conflict(var1, var2)

Describe how two declarations of the same variable disagree. Used to report incompatible duplicates in terms of the attribute which differs rather than by printing both variables and leaving the reader to spot it.

dewpoint_specific_humidity(c, T_dew, p)

Compute near-surface specific humidity [kg/kg] from the dewpoint temperature T_dew [°C] and air pressure p [Pa]. By definition, the actual vapor pressure of the air equals the saturation vapor pressure at the dewpoint, so this delegates to saturation_vapor_pressure — dispatching over ice for T_dew <= 0°C and over liquid water otherwise — and then to vapor_pressure_to_specific_humidity.

diagnose_evaporation_flux(state, model::LandModel)

Return a KernelFunctionOperation that computes the evaporation flux [m/s liquid water] over all grid cells. The kinematic humidity flux is rescaled by the air-to-water density ratio and weighted by the snow-free area fraction.

diagnose_skin_temperature_residual(state, model::LandModel)

Return a KernelFunctionOperation that computes the skin temperature residual over all grid cells.

diagnose_surface_humidity_flux(state, model::LandModel)

Return a KernelFunctionOperation that computes the kinematic surface humidity flux [m/s] over all grid cells.

diffusive_heat_flux(
    i,
    j,
    k,
    grid,
    fields,
    energy::Terrarium.AbstractThermodynamics,
    args...
)

Compute the diffusive (Fourier) heat flux q = -κ ∂T/∂z at the cell face i, j, k, interpolating the medium-specific thermal conductivity to the face.

domains_compatible(d1, d2)

Whether two declarations of the same variable agree on its domain. A declaration without a domain makes no claim and is compatible with any: this is what lets an InputSource, which generally cannot know which domain the variable it feeds belongs to, stay domain-agnostic and adopt whatever the declaring process states. Two different stated domains remain a conflict.

drag_coefficient(i, j, grid, fields, aero::AbstractAerodynamics)

Compute the bulk drag coefficient for heat and moisture transfer at grid cell i, j.

emissivity(i, j, grid, fields, albedo::ConstantAlbedo)

Return the surface emissivity at grid point i, j.

emissivity(i, j, grid, fields, ::AbstractAlbedo)

Return the current emissivity at the given indices.

energy_to_temperature!(
    out,
    i,
    j,
    grid,
    fields,
    ::SnowEnergyTemperatureClosure,
    snow::SingleLayerSnow{NF},
    constants::PhysicalConstants
) where NF

Recover the snow temperature and liquid water fraction from the depth-integrated energy at grid cell i, j.

energy_to_temperature(
    ::FreeWater,
    U::NF,
    ρLθ::NF,
    C::NF
) where NF

Calculate the inverse enthalpy function given the internal energy U (J/m³), volumetric latent heat content ρLθ (J/m³), and volumetric heat capacity C (J/m³/K) under the free water freezing characteristic.

enforce_snow_constraints!(
    state,
    grid,
    snow::Terrarium.AbstractSnow
)

Enforce the constraints on the snow state that are not part of its tendencies: clip negative snow water equivalent to zero and remove the energy of a pack without snow. Called by closure! of models with snow before the snow closure, so that the closure sees an admissible state.

The clipping does not conserve mass; the water it adds to a pack that was overdrawn by an explicit step is not taken from anywhere.

explicit_stage!(state, grid, timestepper, Δt, names)

Apply one explicit stage u += ∂u∂t*Δt to the prognostic fields of state listed in names, using their corresponding tendencies. It is the building block of the explicit timesteppers: ForwardEuler applies one stage per step, and multi-stage schemes such as Heun apply several. Where necessary, the per-field method explicit_stage!(field, tendency, grid, timestepper, Δt) can be specialized for schemes that need a different update.

Accumulate tendency*Δt in the given prognostic field. This method can be overridden by specialized timestepping schemes as needed.

explicit_step!(integrator::ModelIntegrator, timestepper::AbstractTimeStepper, Δt, names::Tuple)

Integrate the prognostic variables in names over a step of size Δt with the explicit timestepper, starting from the current state. The step evaluates the auxiliary variables and tendencies itself (update_state!) and applies one or more explicit_stage!s; multi-stage schemes such as Heun also apply closure! between their stages. The closure at the end of the step and the clock tick are left to the caller: the single-timestepper timestep!, or an AbstractIMEX timestepper, which calls it on its explicit partner and then implicit_step! on its implicit partner before applying the closure once. Explicit timesteppers must implement this method.

explicit_step!(
    integrator::ModelIntegrator,
    timestepper::ForwardEuler,
    Δt,
    names::Tuple
)

Forward Euler step u ← u + Δt ∂u∂t for the prognostic variables in names: compute the auxiliary variables and tendencies at the current state, then apply a single explicit_stage!. The closure is left to the caller (see explicit_step!).

explicit_step!(
    integrator::ModelIntegrator,
    timestepper::Heun,
    Δt,
    names::Tuple
)

Heun step for the prognostic variables in names: a forward Euler predictor with the tendencies at the current state, followed by the closure and a corrector u ← u₀ + Δt (∂u∂t₀ + ∂u∂t₁) / 2 with the tendencies ∂u∂t₁ at the predictor state and the end of the step. The closure after the corrector is left to the caller (see explicit_step!).

explicit_timestepper(imex::AbstractIMEX)

Return the explicit sub-stepper of the AbstractIMEX timestepper imex.

fastiterate(f!::F, iters::NTuple{N,Any}...) where {F,N}

Same as fastmap but simply invokes f! on each argument set without constructing a tuple.

fastmap(f::F, iter::NamedTuple...) where {F}

Same as map for NamedTuples but with guaranteed type stability. fastmap is a @generated function which unrolls calls to f into a loop-free tuple construction expression. All named tuples must have the same keys but in no particular order. The returned NamedTuple

fastmap(f::F, iter::NTuple{N,Any}...) where {F,N}

Same as map for NTuples but with guaranteed type stability. fastmap is a @generated function which unrolls calls to f into a loop-free tuple construction expression.

Compute the empirical field capacity of the soil.

field_indices(
    indices::Tuple{Integer, Integer, Integer}
) -> Tuple{Integer, Integer, Integer}

Pad the grid indices to the three indices required to write to a Field. A 2D (XY) solve passes (i, j), but Oceananigans only defines setindex!(::Field, val, i, j, k) for exactly three indices.

file_extension(::Terrarium.NetCDF) -> String

File extension (including the leading dot) used by files of the given FileFormat.

findfirst_z(i, j, condition_func, z_nodes, field)

2D kernel function that finds the first coordinate in z_nodes where condition_func(field[i, j, k]). This implementation performs a linear scan over the z-axis and thus has time complexity O(N_z).

forcing(
    i,
    j,
    k,
    grid,
    clock,
    fields,
    evapotranspiration::Terrarium.AbstractEvapotranspiration,
    ::Terrarium.AbstractSoilHydrology,
    args...
)

Compute and return the evapotranspiration forcing for soil moisture at the given indices i, j, k. The ET forcing is just the ground_evapotranspiration_flux rescaled by the thickness of layer k.

forcing(i, j, k, grid, clock, fields, forcing::AbstractForcing, target::AbstractProcess, args...)

Return the value computed by the given Oceananigans forcing type, which should be an instance of either DiscreteForcing or ContinuousForcing. Note that target and additional args are only included for interface consistency and are not passed through to forcing.

get_albedo(seb)

Return the albedo parameterization associated with the surface energy balance.

get_artifact(name::String) -> String

Retrieve the path to the artifact with the given name. Throws an AssertionError if no artifact with name exists in Artifacts.toml.

get_asset(asset::Terrarium.AbstractLandAsset) -> String

Download (if necessary) the given asset and return the path to its data file within the installed artifact directory. The artifact is installed via get_artifact and the data file is located by its file_extension. Reading the file is left to the caller; use load_asset to read a variable into a data Field wrapped on the asset's native_grid.

get_biogeochemistry(soil::AbstractSoil)

Return the biogeochemistry scheme associated with soil.

get_cache(cache::AbstractTimeStepperCache, timestepper::AbstractTimeStepper)

Return the working cache for timestepper given the model's state.timestepper_cache. For a single timestepper this is the cache itself; an IMEX cache returns the sub-cache matching the timestepper's class.

get_closure(snow::SingleLayerSnow)

get_closure(
    hydrology::SoilHydrology
) -> Terrarium.AbstractSoilWaterClosure

Return the saturation-pressure closure defined by the given hydrology process, or nothing if not defined for the given configuration.

get_constants(model::AbstractModel)::PhysicalConstants

Return the PhysicalConstants associated with the given model.

get_domain(grid, ::Terrarium.Ground)

Return the spatial discretization of grid for the model domain selected by domain, or nothing if grid does not resolve that domain vertically.

get_energy_balance(soil::AbstractSoil)

Return the energy balance scheme associated with soil.

get_field(
    state,
    ::Union{Terrarium.AbstractVariable{name}, Terrarium.Namespace{name}}
) where name

Retrieves the Field from state matching the name of the given variable.

get_fields(state, ns::Terrarium.Namespace)

Retrieves all Fields declared by the given Namespace from state, where state is assumed to correspond to the (nested) StateVariables of the namespace itself.

get_fields(state, vars)

Retrieves all Fields from state matching the names of the given variables. Any Namespaces in vars are resolved recursively and their fields are merged into the returned NamedTuple keyed by namespace name, with the namespace's own fields collected into a nested NamedTuple.

get_fields(state, components...; except)

Retrieves all non-tendency Fields from state defined on the given components.

get_fields(state, queries::Union{Symbol, Pair}...)

Retrieves fields with names given in queries and returns them in a NamedTuple. Each argument in queries can either be a Symbol corresponding to a field/variable defined in the namespace of state or a Pair{Symbol, Tuple} where the key is the child namespace and the value is a tuple of queries from that namespace.

# initialize model state
state = StateVariables(model)
# get the temperature and saturation_water_ice fields
fields = get_fields(state, :temperature, :saturation_water_ice)
# extract temperature as well as variables from a namespace
nested_fields = get_fields(state, :temperature, :namespace => (:subvar1, :subvar2))

get_grid(model::AbstractModel)::AbstractGrid

Return the spatial grid associated with the given model.

get_ground_heat_flux(seb)

Return the ground_heat_flux process from the surface energy balance.

get_hydraulic_properties(hydrology::AbstractSoilHydrology)

Return the soil hydraulic properties defined by the given soil hydrology configuration.

get_hydraulic_properties(
    hydrology::SoilHydrology
) -> Terrarium.AbstractSoilHydraulics

Return the soil hydraulic properties defined by the given hydrology process.

get_hydrology(soil::AbstractSoil)

Return the hydrology scheme associated with soil.

get_initializer(model::AbstractModel)::AbstractInitializer

Return the initializer associated with the given model.

get_radiative_fluxes(seb)

Return the radiative_fluxes component of the surface energy balance.

get_skin_temperature(seb)

Return the skin_temperature process from the surface energy balance.

get_stratigraphy(soil::AbstractSoil)

Return the stratigraphy parameterization associated with soil.

get_swrc(::AbstractUnsatK)

Return the soil water retention curve associated with the given unsaturated hydraulic conductivity scheme.

get_swrc(hydrology::SoilHydrology)

Return the soil water retention curve from the hydraulic_properties associated with the given SoilHydrology configuration.

get_swrc(hydrology::AbstractSoilHydrology)

Return the soil water retention curve from the hydraulic_properties associated with the given soil hydrology configuration.

get_thermal_properties(energy::AbstractThermodynamics)

Return the thermal properties associated with the given energy balance process.

get_timestepper(model::AbstractModel)::AbstractTimeStepper

Return the timestepper associated with the given model. All AbstractModels are required to define a timestepper field holding an AbstractTimeStepper (e.g. ForwardEuler, Heun, or IMEX).

get_turbulent_fluxes(seb)

Return the turbulent_fluxes component of the surface energy balance.

ground_domain(grid)

Return the spatial discretization of the ground (soil) domain of grid. Grids which are not land grids are their own ground discretization.

ground_evaporation_conductance(
    ::BareGroundEvaporation,
    β,
    rₐ
)

ground_evaporation_conductance(
    ::PALADYNCanopyEvapotranspiration,
    β,
    rₐ,
    rₐ_can
)

Compute the ground evaporation conductance from the given resistance factor β and aerodynamic resistances.

ground_evaporation_resistance_factor(i, j, grid, fields, :AbstractGroundEvaporationResistanceFactor, args...)

Compute the resistance factor against ground evaporation [-] based on the current state and implementation-specific process dependencies in args.

ground_evapotranspiration_flux(i, j, grid, fields, ::AbstractEvapotranspiration)

Return the total ground evapotranspiration flux [m/s], i.e. ground evaporation + plant transpiration, at cell i, j based on the current state.

ground_heat_flux(i, j, grid, fields, ::AbstractGroundHeatFlux)

Return the current ground heat flux at the given indices.

ground_thermal_interface(
    i,
    j,
    grid,
    fields,
    skinT::ImplicitSkinTemperature
) -> Tuple{Any, Any, Any}

Return the conduction target (Tg, κ, Δz) for the snow-free ground: the uppermost ground layer's ground_temperature, the assumed surface conductivity κₛ, and the top ground cell thickness. This is always the unblended ground-only target, used both when there is no snow and (weighted by 1 - f_snow) for the bare-ground share of the skin-temperature solve when there is.

humidity_flux(::Terrarium.AbstractEvapotranspiration, Δq, g)

Compute an evapotranspiration flux (m/s, positive upwards) as the product of a vapor conductance g (m/s) and a specific humidity difference Δq (kg/kg).

Compute (variably saturated) hydraulic conductivity based on the given hydraulic properties, soil water retention curve (SWRC), and volumetric fractions.

hydraulic_conductivity(
    i,
    j,
    k,
    grid,
    fields,
    hydrology,
    strat,
    bgc
)

Compute the hydraulic conductivity at the center of the grid cell i, j, k.

ice_impedance_factor(Ω, f)

Ice impedance factor I_ice = 10^(-Ω (1 - f)) [9], which reduces the hydraulic conductivity of partially frozen soil to account for the blocking of water-filled pores by ice, given the scaling parameter Ω and the unfrozen fraction f of the pore water.

References

• [9] Westermann et al., Geoscientific Model Development (2023)

implicit_closure!(state, model::AbstractModel, ::Iteration{name, target, final})

Process hook of the implicit solver. Re-evaluate the closure of variable name (and the auxiliary variables the next linearization needs) after implicit_update!, for the process that owns the variable only. A model typically forwards it to the regular process closure! of that process; only a closure with once-per-step side effects (e.g. routing excess water to runoff) needs its own method, which must restrict them to final = true (see the soil hydrology for an example).

implicit_closure!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    ::Terrarium.Iteration{name, target, false}
) where {NF, name, target}

Apply the saturation–pressure closure without its once-per-step effects after an intermediate iteration of the implicit step (compute_saturation_closure! without runoff routing or heat transfer), since every iteration restarts from the explicit stage, then update the hydraulic conductivity.

implicit_closure!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    ::Terrarium.Iteration{name, target, true}
) where {NF, name, target}

Apply the regular saturation–pressure closure! after the final iteration of the implicit step, which routes the surface excess to runoff and moves the heat of the redistributed water, then update the hydraulic conductivity.

implicit_operator(
    snow::SingleLayerSnow,
    state,
    grid,
    ::Val{:snow_water_equivalent}
) -> DiagonalOperator

Allocate the DiagonalOperator of the implicit meltwater drainage of snow_water_equivalent, solved in the liquid fraction snow_liquid_fraction (see implicit_variable).

Drainage removes liquid water at the zero-enthalpy reference (liquid water at 0 °C), so it changes neither the snow energy Ū nor the ice content W_i = (1 − θ) W = −Ū / (ρ_w L_f), clamped to [0, W]. With Ū⁺ and W⁺ = Wⁿ + Δt G_exp from the explicit stage (snowfall, rain, and sublimation), the backward Euler step of the drainage M_r(θ, W) (compute_meltwater_outflow) is a scalar equation per column in the liquid fraction θ, via W = W_i / (1 − θ):

G(θ) = W_i / (1 − θ) + Δt M_r(θ, W_i / (1 − θ)) − W⁺ = 0.

Each Picard iteration takes one Newton step: linearizing G about the iterate θᵐ, G(θ) ≈ G(θᵐ) + G′(θᵐ) (θ − θᵐ), and setting this to zero gives the diagonal system G′(θᵐ) θᵐ⁺¹ = G′(θᵐ) θᵐ − G(θᵐ) (see the snow mass balance docs for the convergence argument). Each iteration calls the hooks:

1. linearize!: the diagonal G′(θᵐ) and right-hand side G′(θᵐ) θᵐ − G(θᵐ), at the liquid fraction θᵐ of the current iterate Wᵐ, bounded to [0, θ⁺];
2. solve!: θᵐ⁺¹ in snow_liquid_fraction;
3. implicit_update!: Wᵐ⁺¹ = W_i / (1 − θᵐ⁺¹) with min(L_c, θ⁺) ≤ θᵐ⁺¹ ≤ θ⁺, since the iteration starts from Wⁿ, which may lie below the root, and meltwater_flux ← (W⁺ − Wᵐ⁺¹) / Δt, the drainage passed on to the surface hydrology in the next step;
4. implicit_closure!: the snow closure at the new iterate.

A fully liquid pack (W_i = 0, θ = 1) drains at min(K_sat, (1 − L_c) W / τ_m), so its step has the exact solution max(W⁺ − Δt K_sat, W⁺ / (1 + Δt (1 − L_c) / τ_m)), independent of the iterate.

implicit_operator(
    energy::SoilThermodynamics,
    hydrology::Terrarium.AbstractSoilHydrology,
    state,
    grid
)

Return the implicit operator of the soil internal energy: a TridiagonalOperator on ground_domain(grid) for heat conduction in temperature, with a face velocity for upwind advection of heat by the water flux if the hydrology lets the water flow (SoilHydrology with RichardsEq). Validates on the host that the freeze curve is FreeWater (the only one with an implemented regularized apparent heat capacity) and that temperature has no value boundary conditions (not yet supported). Flux boundary conditions on internal_energy and gradient boundary conditions on temperature are supported; the latter prescribe the diffusive flux -κ ∂T/∂z through the boundary face, which enters the explicit tendency via compute_boundary_energy_tendency.

The operator is the backward Euler step of ∂U/∂t = ∂z(κ ∂z T) - ∂z(C_w T j_w) + G_exp in temperature, with first-order upwind advection by the vertical water flux j_w (T_up the temperature of the upwind cell of a face) and a conservative update of the energy (see the soil energy docs for the derivation). With the Picard solver, for iterations m = 0, …, M-1, starting from U⁰ = Uⁿ (plus any heat moved by the closure of the water step) and T⁰ ← closure(U⁰):

1. linearize!: Tᵐ ← closure(Uᵐ), lagged κᵐ (faces), regularized apparent heat capacity C_appᵐ (centers), advective face velocity C_w j_w (faces), and rhs = C_appᵐ Tᵐ - Uᵐ + U⁺ with U⁺ = U⁰ + Δt G_exp;
2. solve!: C_appᵐ Tᵐ⁺¹ - Δt ∂z(κᵐ ∂z Tᵐ⁺¹) + Δt ∂z(C_w j_w T_upᵐ⁺¹) = rhs;
3. implicit_update!: Uᵐ⁺¹ = U⁺ + Δt [∂z(κᵐ ∂z Tᵐ⁺¹) - ∂z(C_w j_w T_upᵐ⁺¹)] through the interior faces only;
4. implicit_closure!: Tᵐ⁺¹ and the liquid water fraction ← closure(Uᵐ⁺¹).

The solver calls the hooks through the SoilModel/LandModel forwarding methods and the SoilEnergyWaterCarbon methods in soil_coupled.jl.

implicit_operator(model::AbstractModel, state, names::Tuple)

Return a newly allocated AbstractLinearOperator for the implicitly routed prognostic variables names of model (e.g. a SplitOperator of per-variable operators in solve order). Called once on the host at initialize; methods throw an ArgumentError for unsupported configurations, as does the fallback. SoilModel and LandModel forward to implicit_operator(soil::AbstractSoil, state, grid, names).

implicit_operator(
    soil::Terrarium.AbstractSoil,
    state,
    grid,
    names::Tuple{Vararg{Symbol}}
) -> Terrarium.SplitOperator

Return the implicit operator of the soil prognostic variables names on grid: a SplitOperator of per-variable operators that solves saturation_water_ice before internal_energy, whatever the order of names. The heat step then advects heat with the water flux of the current step (water_flux, see advective_heat_flux) and sees the end-of-step saturation, while the Richards step sees the beginning-of-step temperature and liquid water fraction (the effect of frozen ground on the hydraulic conductivity is lagged by one step). Other names follow in their given order and are rejected by the per-variable fallback.

implicit_operator(
    hydrology::Terrarium.AbstractSoilHydrology,
    state,
    grid
)

Fallback for hydrology processes without an implicit operator; throws an ArgumentError.

implicit_operator(
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    state,
    grid
) where NF

Allocate the TridiagonalOperator of saturation_water_ice on ground_domain(grid). Throws an ArgumentError unless the closure is the SoilSaturationPressureClosure and pressure_head has no value boundary conditions.

The operator is the backward Euler step of the Richards equation in the total pressure head Ψ, after the mass-conservative modified Picard iteration of [3]: with porosity φ, saturation ξ, and ξ⁺ = ξⁿ + Δt G_exp from the explicit stage, the water content is linearized about the iterate m, φ ξᵐ⁺¹ ≈ φ ξᵐ + Cᵐ (Ψᵐ⁺¹ - Ψᵐ) with the specific moisture capacity Cᵐ and the hydraulic conductivity Kᵐ lagged, giving the tridiagonal system Cᵐ Ψᵐ⁺¹ - Δt ∂z(Kᵐ ∂z Ψᵐ⁺¹) = Cᵐ Ψᵐ - φ (ξᵐ - ξ⁺) (see the soil hydrology docs for the derivation). Each Picard iteration m = 0, …, M-1, starting from ξ⁰ = ξⁿ, calls the hooks:

1. linearize!: hydraulic conductivity, lagged Kᵐ (faces) and Cᵐ (centers, bounded below by minimum_moisture_capacity so that frozen saturated cells stay regular), and the right-hand side;
2. solve!: the system above for Ψᵐ⁺¹;
3. implicit_update!: ξᵐ⁺¹ = ξ⁺ + Δt ∂z(Kᵐ ∂z Ψᵐ⁺¹) / φ through the interior faces, which conserves water exactly, and water_flux ← -Kᵐ ∂z Ψᵐ⁺¹, the flux with which the heat step advects heat;
4. implicit_closure!: redistribution of over- and undersaturation, water table, Ψᵐ⁺¹, and hydraulic conductivity. Since every iteration restarts from ξ⁺, only the final one routes surface excess water to the runoff pool and moves the heat of redistributed water.

The solver calls the hooks through the SoilModel/LandModel forwarding methods and the SoilEnergyWaterCarbon methods in soil_coupled.jl.

References

• [3] Celia, Bouloutas, and Zarba, Water Resources Research 26, 1483–1496 (1990)

implicit_operator(snow::AbstractSnow, state, grid, ::Val{name})

Return the implicit operator of the snow prognostic variable name: a DiagonalOperator for the meltwater drainage of snow_water_equivalent (see linearize!). The fallback throws an ArgumentError.

implicit_operator(soil::AbstractSoil, state, grid, ::Val{name})

Return the implicit operator of the soil prognostic variable name on grid. SoilEnergyWaterCarbon returns that of the energy process for internal_energy and of the hydrology process for saturation_water_ice; this fallback throws an ArgumentError.

implicit_step!(integrator::ModelIntegrator, timestepper::AbstractTimeStepper, Δt, names::Tuple)

Apply the implicit part of a time step of size Δt with the implicit timestepper to the variables in names. Called by an AbstractIMEX timestepper after the explicit_step! of its explicit partner, so each prognostic field in names already holds u⁺ = uⁿ + Δt G_exp, where Δt G_exp = u⁺ - uⁿ is the explicit increment, which contains only the explicit terms. u⁺ is the right-hand side of the implicit system, which this method solves in place.

This is unrelated to (and does not extend) Oceananigans.TimeSteppers.implicit_step!.

Implicit timesteppers must implement this method; the fallback throws an ArgumentError, and implicit_step!(integrator, ::Nothing, Δt, names) is a no-op.

implicit_step!(
    integrator::ModelIntegrator,
    timestepper::BackwardEuler,
    Δt,
    names::Tuple
)

Apply the implicit Euler step to the prognostic variables in names. On entry, each prognostic field holds u⁺ = uⁿ + Δt G_exp and previous_state of the BackwardEulerCache holds uⁿ (see update_cache!). The step stores Δt G_exp = u⁺ - uⁿ of every variable in the explicit_stage of its ImplicitSolverWorkspace, resets every prognostic field to the initial iterate uⁿ, and then solves all variables with the solver-level implicit_step! (see Picard) using the operator and workspaces of the cache. Since all variables are reset before any is solved, processes see the other implicit variables at uⁿ (or at uⁿ⁺¹ once solved), never at u⁺. On exit, the prognostic fields hold uⁿ⁺¹ and the closure-derived auxiliary variables are consistent with them.

implicit_step!(
    solver::Picard,
    operator::AbstractLinearOperator,
    workspace::Terrarium.ImplicitSolverWorkspace,
    state,
    model::Terrarium.AbstractModel,
    Δt
)

Solve the implicit system of the variable of workspace with the Picard solver and the variable's linear operator. On entry, the prognostic field in state holds the initial iterate u⁰ (uⁿ for BackwardEuler, plus any change δ made by the closures of variables solved before it) and workspace.explicit_stage holds Δt G_exp. The step forms u⁺ = u⁰ + Δt G_exp in workspace.explicit_stage (keeping δ), then for each of the solver.iterations iterations calls linearize! (lagged coefficients of operator and workspace.rhs), solve! (for the target field of the Iteration label workspace.iteration), implicit_update! (conservative update from u⁺), and implicit_closure!, passing the label with its 1-based number and final = true in the last iteration only. On exit, the prognostic field holds uⁿ⁺¹ and the closure variables are consistent with it.

implicit_step!(
    solver::Picard,
    operator::Terrarium.SplitOperator,
    workspaces::NamedTuple,
    state,
    model::Terrarium.AbstractModel,
    Δt
)

Solve the implicit system of a SplitOperator with the Picard solver by operator splitting: the single-variable method below is applied to each component of operator in order, with the ImplicitSolverWorkspace of that variable from workspaces. Each variable is solved at the latest state of the variables solved before it and the initial iterate of those after it, and forms its explicit stage from its current field, as in a sequential split.

implicit_step!(solver::AbstractImplicitSolver, operator::AbstractLinearOperator, workspaces::NamedTuple, state, model::AbstractModel, Δt)
implicit_step!(solver::AbstractImplicitSolver, operator::AbstractLinearOperator, workspace::ImplicitSolverWorkspace, state, model::AbstractModel, Δt)

Solve the implicit system of model with solver and operator (from implicit_operator), using the hooks linearize!, solve!, implicit_update!, and implicit_closure!.

The first form solves all routed variables (workspaces is a NamedTuple of ImplicitSolverWorkspaces keyed by variable name); the second solves the single variable of workspace, e.g. one component of a SplitOperator. On entry, the prognostic fields hold the initial iterate and each explicit_stage holds Δt G_exp; on exit, the prognostic fields hold the solution and their closure variables are consistent with it. The fallbacks throw an ArgumentError.

implicit_timestepper(imex::AbstractIMEX)

Return the implicit sub-stepper of the AbstractIMEX timestepper imex.

implicit_update!(state, model::AbstractModel, operator::AbstractLinearOperator, explicit_stage, Δt, ::Iteration{name})

Process hook of the implicit solver. Update the prognostic variable name in state conservatively from the solved target and explicit_stage, using the fluxes implied by operator.

implicit_update!(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    operator::DiagonalOperator,
    explicit_water,
    Δt,
    ::Terrarium.Iteration{:snow_water_equivalent}
)

Update snow_water_equivalent from the solved liquid fraction and store the drainage of the step in meltwater_flux (compute_implicit_snow_water_update!).

implicit_update!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    operator::TridiagonalOperator,
    explicit_energy,
    Δt
)

Conservative update of the soil internal energy, Uᵐ⁺¹ = U⁺ + Δt [∂z(κ ∂z Tᵐ⁺¹) - ∂z(C_w j_w T_upᵐ⁺¹)], from the implicit diffusive and advective fluxes through the interior faces (compute_implicit_energy_update), with U⁺ = explicit_energy and the lagged face conductivity κ and face velocity C_w j_w of operator.

implicit_update!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    operator::TridiagonalOperator,
    explicit_saturation,
    Δt
) where NF

Update saturation_water_ice conservatively from the solved pressure head (compute_implicit_saturation_update) and store the interior Darcy fluxes in water_flux (compute_implicit_water_flux).

implicit_variable(var::AbstractVariable)

Name of the field that the linear system of the implicitly integrated variable var is solved for if var has no closure: the variable itself by default. Processes may linearize in a derived field instead, e.g. the snow meltwater drainage of snow_water_equivalent in snow_liquid_fraction. Not exported.

Initialize input variables from the given input sources. The scope corresponds to the path of namespace names from the root namespace to state and is used to match namespaced input sources to their target variables; see varpath.

initialize!(state, model::AbstractModel)

Initialize all variables defined in state which are defined by model. This defaults to simply calling initialize!(state, model, get_initializer(model)).

initialize!(state, model::AbstractModel, initializer::AbstractInitializer)

Initialize the model state variables using the corresponding initializer. This method only needs to be implemented if initialization routines are necessary in addition to direct field/variable initializers.

initialize!(state, grid, process::AbstractProcess, args...)

Initialize all state variables associated with the given process on grid. Implementations of AbstractProcess may define additional args that correspond to different process coupling interfaces.

initialize!(inputs, grid, clock, fields, ::InputSource)

Initializes the input source, writing into the input inputs fields at model start. The grid and the full model state fields (read-only) are provided so that sources may compute their inputs from the grid geometry or other state variables. Namespace routing is handled by the caller (see [InputSources]). Default implementation does nothing.

initialize!(state, grid, ::ImplicitSkinTemperature, args...)

Seed the prognostic skin_temperature with the current ground_temperature so the implicit nonlinear solve starts from a physically sensible guess close to the root.

initialize!(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    args...
)

Initialize the snow internal energy by evaluating the inverse closure (temperature → energy). Assumes snow_temperature and snow_water_equivalent have already been initialized.

initialize!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants
)

Initialize the soil energy, water, and carbon state variables on grid given the parameter values in constants.

initialize!(
    state,
    grid,
    hydrology::SoilHydrology,
    soil::Terrarium.AbstractSoil,
    args...
)

initialize!(
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    args...
)

initialize!(state, grid, seb::SurfaceEnergyBalance, args...)

Resets the simulation clock to its initial time (the time of the clock when the integrator was constructed) and calls initialize!(state, model) on the underlying model which should reset all state variables to their values as defined by the model initializer.

initialize!(
    field::Oceananigans.Fields.AbstractField,
    grid::Oceananigans.Grids.AbstractGrid,
    init::Terrarium.AbstractFieldInitializer
) -> Union{Nothing, Field{LX, LY, LZ, O, G, I, D, T} where {LX, LY, LZ, G<:Union{Nothing, Oceananigans.Grids.AbstractGrid}, T, O, I, D}, FieldTimeSeries{LX, LY, LZ, <:Any, <:Oceananigans.OutputReaders.AbstractInMemoryBackend, I, D, G, T} where {LX, LY, LZ, G<:Union{Nothing, Oceananigans.Grids.AbstractGrid}, T, I, D}, Oceananigans.OutputReaders.GPUAdaptedFieldTimeSeries{LX, LY, LZ, <:Any, <:Oceananigans.OutputReaders.AbstractInMemoryBackend, T} where {LX, LY, LZ, T}}

Initialize field on grid with the given field initializer. The default implementation evaluates init(coords...) at the node coordinates of field.

initialize!(
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    args...
) where NF

initialize!(state, inits::NamedTuple{names}) where names

Initialize the state with Field initializers (any valid argument to set!) in inits.

initialize(integrator, params; clock, inputs)

Reconstruct the given integrator using the same underlying model populated with the given params. The clock and inputs can also optionally be updated via their respective keyword arguments.

initialize(
    vars::OrderedCollections.OrderedDict{Symbol, <:Terrarium.AbstractVariable},
    grid::Oceananigans.Grids.AbstractGrid,
    clock::Clock,
    fields::OrderedCollections.OrderedDict{Symbol, Oceananigans.Fields.AbstractField},
    boundary_conditions::NamedTuple
) -> OrderedCollections.OrderedDict{Symbol, Oceananigans.Fields.AbstractField}

Initialize Fields on grid for each of the variables in the given OrderedDict vars. Any predefined boundary_conditions and fields will be passed through to initialize for each variable.

initialize(timestepper::AbstractTimeStepper, state, progvars, model)

Allocate the time stepper cache for timestepper against the given state. progvars is the named tuple of prognostic variable metadata and model the owning AbstractModel; both are needed e.g. by AbstractIMEX timesteppers to resolve each variable's timestepping class.

initialize(
    var::Terrarium.AbstractVariable,
    grid::Oceananigans.Grids.AbstractGrid,
    clock::Clock,
    fields::OrderedCollections.OrderedDict{Symbol, Oceananigans.Fields.AbstractField},
    boundary_conditions::NamedTuple
)

Initialize a Field on grid based on the given var metadata. The named tuple of boundary_conditions should follow the standard convention of (var1 = (; top, bottom, ...), var2 = (; top, bottom, ...)). If user_fields contains a Field matching the name of var, this field will be directly returned. Otherwise check the accumulated dict for fields from previous groups. Otherwise, the new Field is constructed using the given boundary_conditions.

initialize(
    var::Terrarium.AuxiliaryVariable,
    grid::Oceananigans.Grids.AbstractGrid,
    clock::Clock,
    fields::OrderedCollections.OrderedDict{Symbol, Oceananigans.Fields.AbstractField},
    boundary_conditions::NamedTuple
)

Initialize a Field on grid for the given AuxiliaryVariable.

initialize(
    timestepper::BackwardEuler{NF},
    state::StateVariables,
    progvars,
    model
) where NF -> Terrarium.BackwardEulerCache{<:Any, NamedTuple{names, T}, <:AbstractLinearOperator, NamedTuple{names1, T1}} where {names, T<:Tuple, names1, T1<:Tuple}

Allocate the BackwardEulerCache for the prognostic variables in progvars, which under an AbstractIMEX are exactly those routed to the implicit sub-stepper. Per-variable fields and ImplicitSolverWorkspaces are allocated with similar from the prognostic fields of state, the operator with implicit_operator(model, state, names), and the target of each variable with resolve_implicit_target.

Throws an ArgumentError if BackwardEuler is not the implicit sub-stepper of the model's AbstractIMEX timestepper, if any variable in progvars has a non-scalar element type (the implicit solvers are scalar only), if no implicit operator supports the configuration of the variables, or if a target cannot be resolved.

initialize(::AbstractTimeStepper, state)

Initialize and return the AbstractTimeStepperCache holding any intermediate fields/state required by the given timestepper. Time steppers that need no working state fall back to the default implementation, which returns an EmptyCache.

initialize(model::Terrarium.AbstractModel{NF}; ...) where NF
initialize(
    model::Terrarium.AbstractModel{NF},
    params;
    clock::Clock,
    inputs::InputSource,
    boundary_conditions,
    initializers,
    fields
) where NF

Creates and initializes a ModelIntegrator for the given model with input variables populated by the given inputs and optionally params . InputSources can be specified via the inputs keyword argument. This method allocates all necessary Fields for the state variables and subsequently calls initialize!(::ModelIntegrator).

Note that this method is not type stable and thus should not be called from Enzyme autodiff. To reinitialize the model for an existing state, use initialize!(state, model).

See the docstring for initialize(::AbstractModel) for further details.

input(name, loc; default, units, eltype, bounds, desc)

Convenience constructor method for InputVariable.

input_fields(state, components...)

Retrieves all Fields from state corresponding to input variables defined on the given components.

Helper method that selects only input variables declared on obj.

invclosure!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants
)
invclosure!(
    state,
    grid,
    soil::SoilEnergyWaterCarbon,
    constants::PhysicalConstants,
    surface_hydrology::Union{Nothing, Terrarium.AbstractSurfaceHydrology}
)

Compute the inverse closure mapping for soil hydrology and energy, in that order.

See closure! for the role of the optional surface_hydrology process.

invclosure!(
    state,
    grid,
    closure::SnowEnergyTemperatureClosure,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    args...
)

Inverse closure: compute the prognostic snow_energy from a prescribed snow_temperature (used for initialization). snow_water_equivalent must be initialized first, since the snow depth d_snow enters the depth-integrated energy.

invclosure!(state, grid, closure::AbstractClosureRelation, process::AbstractProcess, args...)

Apply the inverse of closure for the process proc on the given grid with additional implementation-specific args defined by the coupling interface for the process type.

invclosure!(state, grid, proc::AbstractProcess, args...)

Apply the inverse closure mappings for the process proc on the given grid with additional implementation-specific args defined by the coupling interface for the process type. The default implementation calls invclosure! with args for each closure returned by closures.

invclosure!(state, model::AbstractModel)

Apply the inverse of all closure relations defined for the given model.

invclosure!(
    state,
    grid,
    closure::SoilSaturationPressureClosure,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil
) where NF
invclosure!(
    state,
    grid,
    closure::SoilSaturationPressureClosure,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    runoff::Union{Nothing, Terrarium.AbstractSurfaceRunoff},
    args...
) where NF

Computes saturation_water_ice from the current pressure_head state.

See closure! for the role of the optional runoff process.

iterations(::Terrarium.NewtonSolver{NF, N}) where {NF, N}

Number of Newton iterations performed by the given solver.

kelvin_to_celsius(c::ThermodynamicConstants, T)

Convert the given temperature in °C to Kelvin based on the constant temperature_water_freeze.

kernel(func, args...; clock = false)

Return a KernelFunction that lazily constructs a KernelFunctionOperation from the given func and tuple of args when invoked, i.e:

ctor = kerenl(my_function, arg1, arg2)
...
kfo = ctor(grid, clock, fields)

This is intended to be used as a constructor for AuxiliaryVariables:

myvar(i, j, k, grid, fields) = clamp(fields.x[i, j, k], zero(eltype(grid)), one(eltype(grid)))

auxvar = auxiliary(:myvar, Ground(XYZ()), kernel(myvar))

kernel_operation(
    func,
    state,
    grid,
    args...;
    location,
    with_clock
)

Return a KernelFunctionOperation over grid with the given state and trailing args.

latent_heat_flux(i, j, grid, fields, ::AbstractTurbulentFluxes)

Return the current latent heat flux at the given indices.

launch_kernel!(
    grid,
    workspec,
    kernel!,
    first_arg,
    other_args...;
    kwargs...
)

Launch kernel! over grid with the Oceananigans workspec. grid is passed to the kernel as its second argument, after first_arg.

linearize!(
    operator::DiagonalOperator,
    rhs,
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    explicit_water,
    Δt,
    ::Terrarium.Iteration{:snow_water_equivalent}
)

Fill the diagonal and right-hand side of the Newton step for the liquid fraction of the meltwater drainage (compute_implicit_meltwater_coefficients!), given the explicit stage W⁺ = explicit_water.

linearize!(
    operator::TridiagonalOperator,
    rhs,
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    explicit_energy,
    Δt
)

Linearize the implicit heat operator of the soil internal energy at the current iterate (linearly implicit backward Euler in temperature, conservative update of the internal energy). On entry, explicit_energy holds U⁺ = U⁰ + Δt G_exp, where G_exp contains the flux boundary conditions (e.g. SoilHeatFlux and GeothermalHeatFlux), the diffusive fluxes implied by gradient boundary conditions on temperature, the advective heat fluxes through the boundary faces, and the heat carried by water crossing the boundaries (compute_boundary_advected_heat!), but not the fluxes through the interior faces. This hook

1. re-evaluates the energy closure, so that temperature and liquid water fraction are consistent with the internal energy and the current saturation (changed by the water step since tⁿ),
2. evaluates the face conductivity κ (compute_implicit_conductivity!), the regularized apparent heat capacity C_app (compute_implicit_capacity!), and, with water flow, the advective face velocity C_w j_w from the water_flux of the step (compute_implicit_advective_velocity!), storing them in the face_conductivity, diagonal_coefficient, and face_velocity of operator,
3. and assembles the right-hand side C_app Tᵐ - Uᵐ + U⁺ of C_app Tᵐ⁺¹ - Δt ∂z(κ ∂z Tᵐ⁺¹) + Δt ∂z(C_w j_w T_upᵐ⁺¹) = C_app Tᵐ - Uᵐ + U⁺ (compute_implicit_energy_rhs).

implicit_update! then updates the internal energy conservatively: for any number of iterations, the change of column-integrated energy is exactly Δt times the net boundary flux in G_exp.

linearize!(
    operator::TridiagonalOperator,
    rhs,
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    explicit_saturation,
    Δt
) where NF

Update the hydraulic conductivity (compute_hydraulics!), fill the lagged capacity and face conductivity of operator (compute_implicit_water_flow_coefficients!), and assemble the right-hand side C Ψᵐ - φ (ξᵐ - ξ⁺) in rhs (compute_implicit_water_flow_rhs), with ξ⁺ = explicit_saturation. [3] recommend at least two iterations.

References

• [3] Celia et al., Water Resources Research 26, 1483–1496 (1990)

linearize!(operator::AbstractLinearOperator, rhs, state, model::AbstractModel, explicit_stage, Δt, ::Iteration{name})

Process hook of the implicit solver. Fill the lagged coefficients of operator and the right-hand side rhs for variable name at the current iterate in state, given the explicit stage explicit_stage (u⁺ = uⁿ + Δt G_exp). The fallback throws an ArgumentError.

liquid_water_fraction(i, j, k, grid, fields, ::AbstractSoilHydrology)

Compute or retrieve the current fraction of unfrozen water in the pore space.

liquid_water_fraction(::FreeWater, U::NF, ρLθ::NF) where NF

Calculate the unfrozen water content from the given internal energy U (J/m³) and volumetric latent heat content ρLθ (J/m³).

litterfall_bc(i, j, grid, clock, fields)

Discrete-form boundary-condition function (see LitterfallFlux), computes incoming litterfall carbon inputs to the soil. Note that module computes infiltration as positive downward, so it is negated here since fluxes are by convention positive upward.

load_asset(
    asset::Terrarium.AbstractLandAsset,
    name::String;
    NF,
    fill_value
)

Download (if necessary) the given asset via get_asset and read the variable name from its data file, returning a suitable Array-type based on the asset's file format. The asset's indices are applied when reading (e.g. to select a single time record), and fill_value replaces missing data; it defaults to NF(NaN).

The underlying read is dispatched to an I/O extension based on the asset's format; load Rasters.jl and NCDatasets.jl to enable reading NetCDF and other raster files.