Private API

ConstitutiveParameters(FT = Oceananigans.defaults.FloatType;
                       gas_constant       = 8.3144598,
                       dry_air_molar_mass = 0.02897,
                       water_molar_mass   = 0.018015)

Construct a set of parameters that define the density of moist air,

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

where $p$ is pressure, $T$ is temperature, $q$ defines the partition of total mass into vapor, liquid, and ice mass fractions, and $Rᵐ$ is the effective specific gas constant for the mixture.

HeatCapacityParameters(FT = Oceananigans.defaults.FloatType;
                       dry_air_adiabatic_exponent = 2/7,
                       water_vapor_heat_capacity = 1859,
                       liquid_water_heat_capacity = 4181,
                       water_ice_heat_capacity = 2100)

Isobaric heat capacities.

_hydrostatic_pressure_from_surface!(
    dev
) -> Union{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.DynamicSize, typeof(NumericalEarth.Atmospheres.cpu__hydrostatic_pressure_from_surface!)}, KernelAbstractions.Kernel{Backend, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.DynamicSize, typeof(NumericalEarth.Atmospheres.gpu__hydrostatic_pressure_from_surface!)} where Backend<:KernelAbstractions.GPU}

Hydrostatically-balanced pressure on the grid of temperature, integrating d(ln p)/dz = −g/(Rᵐ T) upward from surface_pressure (the surface pressure at the orography height) with the moist mixture gas constant Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ, qᵈ = 1 − qᵛ − qᶜ − qⁱ.

Each column is anchored at the grid's bottom face — the terrain surface on a terrain-following grid, the domain bottom on a regular height grid — and surface_pressure is reduced hydrostatically from orography to that face. Anchoring on the surface pressure (rather than interpolating pressure to the node heights) avoids the spurious near-surface density a sub-surface-clamped pressure interpolation produces over high terrain.

surface_pressure and orography are (Nx, Ny) arrays. The moisture fields qᵛ, qᶜ, qⁱ are optional Fields; omitting them gives the dry result (Rᵐ = Rᵈ). The calculation launches one kernel thread per horizontal column, integrating upward through k.

using Oceananigans
using NumericalEarth

grid = RectilinearGrid(size = (1, 1, 4), x = (0, 1), y = (0, 1), z = (0, 1000),
                       topology = (Periodic, Periodic, Bounded))
T = CenterField(grid)
set!(T, 250)

p = hydrostatic_pressure_from_surface(T, fill(1.0e5, 1, 1), fill(0.0, 1, 1);
                                      dry_gas_constant = 287.0,
                                      vapor_gas_constant = 461.0,
                                      gravitational_acceleration = 9.81)
round(Int, interior(p)[1, 1, 1] / 100)  # surface pressure (hPa) at the lowest level

# output
983

The standard unit of atmospheric pressure; 1 standard atmosphere (atm) = 101,325 Pascals (Pa) in SI units. This is approximately equal to the mean sea-level atmospheric pressure on Earth.

Interpolate the atmospheric state onto the ocean / sea-ice grid.

set!(atmosphere::PrescribedAtmosphere; u=nothing, v=nothing, T=nothing, q=nothing, p=nothing)

Set the prescribed atmospheric velocities (u, v), air temperature (T), specific humidity (q), and surface pressure (p), then refresh the interpolated state. Omitted keywords are left untouched. A Number sets a constant in space and time.

read_2d_nemo_variable(ds, name)

Read a 2D variable from a NEMO NetCDF dataset, handling varying dimension layouts: (x, y), (x, y, z), or (x, y, z, t).

read_orca_staggered_mesh(ds, overlap)

Read ORCA horizontal coordinates and metrics.

Supports:

• full NEMO staggered mesh variables (glamt/gphit/e1u/...), and
• approximate reconstruction from T/F coordinates only (glamt/gphit/glamf/gphif) using Tripolar-style spherical metric assumptions.

remove_minor_basins!(z_data, keep_major_basins)

Remove independent basins from the bathymetry data stored in z_data by identifying connected regions below sea level. Basins are removed from smallest to largest until only keep_major_basins remain.

Arguments

• z_data: A 2D array representing the bathymetry data.
• keep_major_basins: The maximum number of connected regions to keep. If Inf is provided then all connected regions are kept.

variable_glossary :: Dict{Symbol, Symbol}

Global map from verbose dataset variable names (the symbols a user passes to Metadata and MetadataSet) to the short model field-name symbols established in docs/src/appendix/notation.md. Used by set!(model, mset) to auto-route mset.eastward_velocity → u = ... etc.

Verbose names absent from this map are silently ignored by set!(model, mset) (they remain fetchable via download(mset) and accessible via mset.<name>). Synonyms across dataset modules (e.g. :u_velocity, :eastward_velocity, :eastward_wind all → :u) are intentional: they serve as domain disambiguators when a single dataset (e.g. ECCO4Monthly) carries both ocean and atmosphere fields.

Every value here is documented in docs/src/appendix/notation.md (or in Breeze.jl's notation for the microphysics symbols).

AbstractStaticBathymetry <: AbstractStaticDataset

Supertype for static, two-dimensional bathymetry datasets (e.g. ETOPO, GEBCO, IBCSO, IBCAO). Adds defaults for the degenerate vertical axis and a variable-agnostic Base.size.

AbstractStaticDataset

Supertype for datasets without a time dimension. Provides default no-op implementations for the date-related interface (all_dates, first_date, last_date).

ColumnInfo{F, I}

Resolved location of a Column extraction inside the file grid. Built once per set_region_data! call by region_info(::Column, …) and captured into _set_region_kernel! as a stack-friendly struct.

• i⁻, i⁺: bracketing longitude indices (i⁺ wraps to 1 across the periodic seam).
• j⁻, j⁺: bracketing latitude indices.
• wx, wy: bilinear blend weights in [0, 1] (0 → at i⁻/j⁻, 1 → at i⁺/j⁺).
• ℑ: interpolation kind, Linear() or Nearest().

DatasetBackend{N, C, I, M} <: AbstractInMemoryBackend{Int}

In-memory backend for a FieldTimeSeries backed by a dataset whose metadata maps each in-memory time index to a file (or subset of a file) on disk. The backend carries

• start, length: sliding-window extents into the metadata
• inpainting: inpainting algorithm used when reading per-file datasets (e.g. NearestNeighborInpainting); nothing for datasets whose native NetCDF already covers the whole target grid (e.g. JRA55).
• metadata: the dataset metadata — its dataset type parameterizes set! dispatch, so per-file and chunked (multi-time-per-file) datasets can coexist under the same backend.

Type parameters N (on-native-grid) and C (cache-inpainted-data) are flags hoisted into the type so that dispatch and Adapt.adapt_structure can act on them without allocating.

DatasetBackend(length, metadata;
               on_native_grid = false,
               cache_inpainted_data = false,
               inpainting = NearestNeighborInpainting(Inf))
DatasetBackend(start, length, metadata; ...)

Construct a DatasetBackend holding length in-memory time indices starting at start (default 1). inpainting = nothing selects the dispatch path for datasets whose files hold multiple time instances (e.g. chunked NetCDF like JRA55), where per-file inpainting is not applicable.

DownloadProgress(total, now; filename="")

NearestNeighborInpainting{M}

A structure representing the nearest neighbor inpainting algorithm, where a missing value is substituted with the average of the surrounding valid values. This process is repeated a maximum of maxiter times, until the field is completely inpainted, or until a pass fills nothing.

Field(metadata::Metadatum, arch=CPU();
      inpainting = default_inpainting(metadata),
      mask = nothing,
      halo = (3, 3, 3),
      cache_inpainted_data = true)

Return a Field on architecture described by metadata with halo size. If not nothing, the inpainting method is used to fill the cells within the specified mask. mask is set to compute_mask for non-nothing inpainting. Keyword argument cache_inpainted_data dictates whether the inpainted data is cached to avoid recomputing it; default: true.

Field(mset::MetadataSet, arch=CPU(); kw...)

Build a NamedTuple of Fields — one per variable in mset, keyed by the verbose dataset variable name. Each value is Field(mset[name], arch; kw...).

Requires mset to hold scalar dates so each mset[name] is a Metadatum; for multi-date sets, build a NamedTuple of FieldTimeSeries per variable, e.g. NamedTuple(name => FieldTimeSeries(mset[name], grid) for name in mset.names).

Field(metadata::Metadatum, grid::AbstractGrid; cache = false, overwrite_cache = false, kw...)

Load metadata on its native grid and interpolate onto grid — the Field analog of FieldTimeSeries(metadata, grid). Keyword arguments are forwarded to the native-grid Field(metadata, arch; …) (e.g. inpainting, mask, halo, cache_inpainted_data).

With cache = true the regridded result is cached to disk and reused by later reads with the same dataset, variable, date, region, target-grid geometry, and read keywords — skipping the native materialization and regrid entirely; with cache = false (default) the cache is disabled entirely and nothing is read or written. The key carries a size/mtime stamp of the local dataset file where one exists, so a re-download invalidates the cache. For streaming datasets with no local file, pass overwrite_cache = true after replacing data upstream: it skips the lookup and overwrites the entry with a freshly regridded result.

FieldTimeSeries(metadata::Metadata [, arch_or_grid=CPU() ];
                time_indices_in_memory = 2,
                time_indexing = Cyclical(),
                inpainting = nothing,
                cache_inpainted_data = true)

Create a FieldTimeSeries from a dataset that corresponds to metadata.

Arguments

• metadata: Metadata containing information about the dataset.

• arch_or_grid: Either a grid to interpolate the data to, or an architecture to use for the native grid. Default: CPU().

Keyword Arguments

• time_indices_in_memory: The number of time indices to keep in memory. Default: 2.

• time_indexing: The time indexing scheme to use. Default: Cyclical().

• inpainting: The inpainting algorithm to use for the interpolation. The only option is NearestNeighborInpainting(maxiter), where an average of the valid surrounding values is used maxiter times.

• cache_inpainted_data: If true, the data is cached to disk after inpainting for later retrieving. Default: true.

download(mset::MetadataSet; kwargs...)

Download every variable in mset. The default is a per-variable loop calling download(mset[name]; kwargs...); backends that support batched multi-variable requests (e.g. the ERA5 pressure-level CDS path) override this to route through a single batched call.

Returns a NamedTuple keyed by the set's variable names, whose values are the results of each per-variable download call (typically the file path(s)).

available_variables(metadata)

Return the available variables in the dataset.

averaging_window(metadatum)

The (start, stop) dates of the interval of time the value in metadatum represents, its time_bnds in CF terms. Defaults to a zero-width window at the stamp itself, which is a value at an instant; a product whose files hold time averages extends this with the interval each average is taken over. Whether the dates label the start of that interval, its end, or its middle varies from product to product. See calendar_month_window for the calendar-month case, whose length varies with the month.

bounding_box_intersects(bounds, bbox::BoundingBox)

Whether bounds, anything with west, south, east, and north fields, overlaps bbox. Used to select the files of a tiled dataset that cover a region.

Both boxes must label longitudes in the same convention, with west < east; a box folded across the antimeridian is not supported.

bounding_box_suffix(region)

Filename suffix identifying the window a regionally-downloaded dataset was cached for: "global" for nothing, and "lon_<west>_<east>_lat_<south>_<north>" for a BoundingBox, so windows of one dataset never collide on disk. Datasets whose products follow a different filename convention define their own suffix instead.

calendar_month_window(metadatum)

The averaging_window of a product whose files hold calendar-month means: the month containing the date metadatum is stamped with. Independent of where in the month the stamp falls, so it also covers products stamped at, say, noon on the first.

centers_to_interfaces(z_centers)

Compute $z$-interfaces (cell faces) from cell center positions. z_centers should be sorted most negative first (deepest first). The top face is placed at 0.0 (sea surface). Interior faces are midpoints between adjacent centers. The bottom face is extrapolated.

Note: the grid's cell centers (midpoints of faces) will approximately but not exactly match the input centers when spacing is irregular.

cmr_granules(short_name, version, bbox; extension = "h5", page_size = 2000, attempts = 3)

Return the download URLs of the short_name / version granules whose footprints intersect bbox, querying NASA CMR page by page until a short page signals the last one. Only URLs ending in extension are collected, and one URL is kept per granule name (preferring a protected data-host endpoint).

cmr_granules_url(short_name, version, bbox; page_size = 2000, page_num = 1)

Build the NASA CMR granule-search URL for page page_num of the product short_name / version whose granules intersect the bbox BoundingBox (encoded W,S,E,N, longitudes in [-180, 180]). CMR search is anonymous; only the granule download itself needs Earthdata credentials.

compute_native_date_range(native_dates, start_date, end_date)

Compute the range of native_dates that fall within the specified start_date and end_date.

compute_mask(metadata::Metadatum, dataset_field,
             mask_value = default_mask_value(metadata),
             minimum_value = -1f5,
             maximum_value = 1f5)

A boolean field where true represents a missing value in the dataset_field.

continue_downwards!(field, mask)

Continue downwards a field with missing values within mask. Cells where mask[i, j, k] == false will be preserved.

conversion_units(metadatum)

Return the units of the source variable in the given dataset referenced by metadatum. These units will be used to apply automatic conversions to standard units for NumericalEarth.

dataset_location(dataset, variable_name)

Return the native field location (LX, LY, LZ) for variable_name in dataset. Defaults to (Center, Center, Center). Only datasets with staggered variables (e.g., ECCO velocity fields) need to extend this.

dataset_variable_name(metadata)

Return the name used for the variable metadata.name in its raw dataset file.

default_download_directory(dataset)

Return the default directory to which dataset is downloaded.

default_horizontal_padding(dataset)

Return the default horizontal padding (degrees) added around a bounding box requested from dataset, providing margin for interpolation stencils at the boundary.

default_region(dataset)

Return the default region used when constructing Metadata/Metadatum for dataset. Defaults to nothing (full domain). Single-location datasets (e.g. a moored buoy) can override this to embed their fixed Column so callers need not repeat the coordinates.

download_cache(key)

Return the directory used to cache key data downloaded by NumericalEarth.

By default this is a Scratch.jl space managed by Julia under the active depot. If the environment variable NUMERICALEARTH_DATA_DIRECTORY is set, data is instead cached under joinpath(ENV["NUMERICALEARTH_DATA_DIRECTORY"], key). This is useful on systems where the Julia depot lives on a small or quota-limited filesystem (e.g. $HOME on HPC clusters), or to share a single cache of large datasets across depots and users.

The variable is read when NumericalEarth is loaded, so it must be set before using NumericalEarth.

download_with_retries(
    url,
    path;
    attempts,
    description,
    kw...
) -> Union{Nothing, AbstractString}

Download url to path, retrying up to attempts times on failure. Each transfer lands beside the destination and is renamed into place once it has arrived in full, so path never holds a partial file.

Extra keyword arguments are forwarded to Downloads.download (progress, downloader, ...).

earthdata_download(url, path; attempts = 3)

Download the NASA Earthdata granule at url to path, authenticating with the EARTHDATA_USERNAME / EARTHDATA_PASSWORD environment variables (register free at https://urs.earthdata.nasa.gov).

earthdata_download_cached(url, cache_dir; attempts = 3)

Download the granule at url into cache_dir unless it is already there, and return its path. Keyed on the granule name, so overlapping regions and sibling variables reuse a granule instead of re-downloading it.

expand_dates(dataset, variable_name, dates)

Return dates, expanding a (start_date, end_date) tuple to the dataset's native dates in that window — the cadence is the dataset's own.

figshare_article_url(id)

Return the figshare API URL for article id, which serves the article's metadata including one entry per file, each with a name and a download_url.

fill_gaps!(fts::FieldTimeSeries; max_gap=6)
fill_gaps!(data::AbstractArray; max_gap=6)
fill_gaps!(data::AbstractVector; max_gap=6)

Fill NaN gaps along the time dimension using linear interpolation, except at the ends of the record, where a gap takes the nearest valid value. For an AbstractArray, the last dimension is assumed to be time, and each spatial column is filled independently. For a FieldTimeSeries, interior(fts) is copied to the CPU, filled in place, and copied back.

Interior gaps longer than max_gap points are left as NaN with a warning.

inpaint_mask!(field, mask; inpainting=NearestNeighborInpainting(Inf))

Inpaint field within mask, using values outside mask. In other words, cells where mask[i, j, k] == true are inpainted and cells where mask[i, j, k] == false are preserved.

Arguments

• field: Field to be inpainted.

• mask: Boolean-valued Field, values where mask[i, j, k] == true are inpainted. nothing leaves field untouched.

• inpainting: The inpainting algorithm to use. The only option is NearestNeighborInpainting(maxiter), which substitutes an average of the valid surrounding values at most maxiter times; a bare Int is taken as that maxiter. Default: NearestNeighborInpainting(Inf).

inpainted_metadata_filename(metadata::Metadatum) -> String

Name of the file holding the inpainted version of metadata. Defaults to the metadatum's own filename with its extension replaced by _inpainted.jld2; datasets whose inpainted products are named differently override this.

inpainted_metadata_path(metadata::Metadatum) -> String

Path to the inpainted version of metadata, alongside the raw file in metadata.dir.

interpolate_physical!(to_field, from_field)

Interpolate from_field onto to_field. Identical to Oceananigans' interpolate! for ordinary grids, but on a MutableVerticalDiscretization (terrain-following) target it samples the source at the target's deformed znode rather than the reference rnode. interpolate! builds its target node from _node, whose vertical component is rnode, so it would place the source at the LAM's reference heights and ignore the terrain — putting the lowest cells below the (clipped) surface of a PressureLevelGrid source.

matching_single_level_dataset(dataset)

Return the single-level (surface) dataset sharing dataset's product family and temporal cadence — e.g. the surface companion of a pressure-level reanalysis product.

metadata_path(mset::MetadataSet)

Return a NamedTuple keyed by the set's variable names whose values are the file paths of each variable's Metadata (a String for single-date sets, a Vector{String} for multi-date sets — matching metadata_path(::Metadata)).

metadata_url(metadata)

Return the URL for the dataset described by metadata. Extended by each dataset module.

missing_value(metadatum)

Return the value used by the underlying dataset to represent missing data. Defaults to missing.

native_region_grid(region::BoundingBox, Δλ, Δφ; pad = 2)

Regular lat/lon raster of cell steps Δλ/Δφ (degrees) covering region, snapped to the global lattice anchored at (-180, -90) and padded by pad cells on each side. Returns (; west, south, Δλ, Δφ, Nx, Ny).

Datasets distributed as vector or tiled files lay out the raster they burn onto with this, so the result is a sub-window of the global lattice Field(::Metadatum) assumes.

native_times(metadata; start_time=first(metadata).dates)

Extract the time values from the given metadata, calculate the time difference from the start_time, and return an array of time differences in seconds.

Argument

• metadata: The metadata containing the date information.

Keyword Argument

• start_time: The start time for calculating the time difference. Defaults to the first date in the metadata.

Times are measured to each sample's window_center, which for a time average is the midpoint of its averaging_window rather than the date its file is stamped with.

netrc_downloader(username, password, machine, dir; verify_ssl = true)

Create a downloader that uses a netrc file to authenticate with the given machine. Pass verify_ssl = false only for hosts serving an untrusted CA certificate, as ecco.jpl.nasa.gov does. This downloader writes the username and password in a file named auth.netrc (for Unix) and auth_netrc (for Windows), located in the directory dir. To avoid leaving the password on disk after the downloader has been used, it is recommended to initialize the downloader in a temporary directory, which will be removed after the download is complete.

For example:

mktempdir(dir) do tmp
    dowloader = netrc_downloader(username, password, machine, tmp)
    Downloads.download(fileurl, filepath; downloader)
end

prescribed_radiation(
    dataset;
    ...
) -> PrescribedRadiation{_A, Float64, _B, SW, LW, S} where {_A, _B, SW<:(FieldTimeSeries{LX, LY, LZ, TI, K, I, D, G, ET, B, χ, Nothing, Nothing, @NamedTuple{}} where {LX, LY, LZ, TI, K, I, D, G, ET, B, χ}), LW<:(FieldTimeSeries{LX, LY, LZ, TI, K, I, D, G, ET, B, χ, Nothing, Nothing, @NamedTuple{}} where {LX, LY, LZ, TI, K, I, D, G, ET, B, χ}), S<:NamedTuple}
prescribed_radiation(
    dataset,
    architecture;
    start_date,
    end_date,
    dir,
    time_indices_in_memory,
    time_indexing,
    ocean_surface,
    sea_ice_surface,
    snow_surface,
    land_surface,
    stefan_boltzmann_constant,
    region,
    other_kw...
) -> PrescribedRadiation{_A, Float64, _B, SW, LW, S} where {_A, _B, SW<:(FieldTimeSeries{LX, LY, LZ, TI, K, I, D, G, ET, B, χ, Nothing, Nothing, @NamedTuple{}} where {LX, LY, LZ, TI, K, I, D, G, ET, B, χ}), LW<:(FieldTimeSeries{LX, LY, LZ, TI, K, I, D, G, ET, B, χ, Nothing, Nothing, @NamedTuple{}} where {LX, LY, LZ, TI, K, I, D, G, ET, B, χ}), S<:NamedTuple}

Return a PrescribedRadiation backed by dataset's downwelling shortwave and longwave FieldTimeSeries.

region (a BoundingBox) restricts the download and the native grid to a sub-domain; the coupled model interpolates the native-resolution radiation onto the exchange grid. Surface radiative properties (albedo, emissivity) default to standard ocean/sea-ice values; pass *_surface = nothing to omit a surface or supply your own SurfaceRadiationProperties (e.g. land_surface for land-only runs).

propagate_horizontally!(inpainting, field, mask [, substituting_field=deepcopy(field)])

Horizontally propagate the values of field into the mask. In other words, cells where mask[i, j, k] == false are preserved, and cells where mask[i, j, k] == true are painted over.

The first argument inpainting is the inpainting algorithm to use in the _propagate_field! step.

propagate_vertically!(field; maxiter = size(field, 3))

Fill the NaN cells of field by repeatedly averaging their valid vertical neighbours, the $z$-axis analogue of propagate_horizontally!.

regrid_from_metadata!(target, metadata; tile_bytes = default_tile_bytes, kw...)

Fill target from metadata, regridding in tiles of at most tile_bytes of source each where the dataset reads windows, and in a single pass otherwise. Remaining keyword arguments go to Field(metadata, arch; …).

regrid_in_tiles!(target, metadata, native, tile_bytes)

Fill target from metadata in pieces holding at most tile_bytes of source each, reading only the window each piece needs.

retrieve_data(metadata)

Retrieve data from netcdf file according to metadata.

retrieve_window(metadata, longitude_indices, latitude_indices)

Return (data, λ, φ) for one horizontal window of metadata's file: the data, always three-dimensional and ascending in latitude and height, together with the cell centers the window spans. Indices are ranges, or Colon to span the whole file axis. Serves the hyperslab of a plain lon/lat NetCDF, north-first files included; datasets stored otherwise extend this.

set_region_data!(target, data, λc, φc, metadata)

Fill the region of target (Field or FieldTimeSeries) implied by metadata.region from data, applying mangling, NaN conversion, and unit conversion in a single GPU-friendly kernel pass.

uncovered_time_gaps(metadata)

The (head, tail) durations in seconds that the averaging_windows of metadata span but its nodes do not: from the start of the first window to the first window_center, and from the last center to the end of the last window. A FieldTimeSeries extrapolates over both. Zero for a product of values at an instant, whose nodes span exactly its samples.

unzip(file) -> String
unzip(file, exdir) -> AbstractString

Extract the zip archive file into the directory exdir, which is created if missing.

validate_time_coverage(metadata, time_indexing)

Check that time_indexing is a valid operation outside the nodes of a window-averaged series.

Nodes sit at window_centers, so the first and last half window of the span the data covers falls outside them — see uncovered_time_gaps. Clamp holds the end values there; Cyclical fills them from the far end of the record; Linear holds the first value before the first node but extrapolates without bound past the last, so it can leave the range of the data.

A single sample has no far end and nothing to interpolate between (it is constant in time) so it is left alone whatever the scheme.

window_center(metadatum)

The date at the middle of the averaging_window of metadatum, where a time average equals the value of the field itself and so where a linearly interpolating FieldTimeSeries has to place it. The stamp itself for a value at an instant.

window_span(metadata)

The span in seconds from the start of the first averaging_window in metadata to the end of the last, which is the period over which a window-averaged series repeats: its windows tile that span without gaps or overlaps. nothing for a product of values at an instant, whose samples are points and so leave the period to be inferred from the node spacing.

windowed_retrieval(dataset)

Whether dataset serves a horizontal subwindow of its file through retrieve_window, which lets Field(metadatum, grid) regrid it in tiles. A dataset that overrides interpolate_physical! regrids by its own scheme and leaves this false. Default false.

z_interfaces(dataset)

Return an array with the vertical interfaces ($z$-faces) of the dataset that metadata corresponds to.

set!(model, mset::MetadataSet, names=keys(variable_glossary))

Route variables from mset to set!(model; kwargs...), translating verbose dataset names to short model field-names via variable_glossary. Only the intersection of names and mset.names is forwarded.

The default names is every glossary key — fine for permissive models. Models that throw on unknown kwargs (HydrostaticFreeSurfaceModel, SeaIceModel) override the 2-argument form to pass a narrower names, letting a single multi-component MetadataSet drive both an ocean and a sea-ice model. Each model's override of set!(model, ::MetadataSet) filters the same mset down to the variables that model knows how to consume:

using NumericalEarth
using Oceananigans
using Statistics
using Dates

grid = LatitudeLongitudeGrid(size = (60, 30, 5),
                             longitude = (-180, 180),
                             latitude  = (-60, 60),
                             z = (-5000, 0),
                             halo = (7, 7, 7))

ocean = ocean_simulation(grid)

mset = MetadataSet(:temperature, :salinity;
                   dataset = ECCO4Monthly(),
                   date    = DateTime(1993, 1, 1))

# Ocean override routes :temperature → T, :salinity → S; sea-ice vars are
# filtered out. A `set!(sea_ice.model, mset)` call against the same `mset`
# would route :sea_ice_thickness → h, :sea_ice_concentration → ℵ.
set!(ocean.model, mset)

T = ocean.model.tracers.T
(min = round(minimum(T), digits = 2),
 max = round(maximum(T), digits = 2),
 mean = round(mean(T), digits = 2))

# output
(min = -1.06, max = 21.41, mean = 3.3)

set!(fields::NamedTuple, mset::MetadataSet)

Set each fields[name] from the corresponding mset[name]. The NamedTuple's keys must be a subset of the set's variable names; extra fields are ignored.

This is the explicit form that takes verbose dataset names on both sides — no glossary translation. For the auto-routing form (short model field-names), see set!(model, ::MetadataSet).

asterged_decode_emissivity(DN)

Decode a raw /Emissivity/Mean digital number to a Float32 emissivity: fill value −9999 maps to NaN, otherwise ε = 0.001 · DN clamped to [0, 1]. The clamp is needed because ASTER's Temperature/Emissivity Separation retrieval carries no ε ≤ 1 constraint, and the product stores band emissivities above unity.

asterged_decode_uncertainty(DN)

Decode a raw /Emissivity/SDev digital number to a Float32 emissivity standard deviation: fill value −9999 maps to NaN, otherwise σ = 1e-4 · DN. Note the scale is 10× smaller than asterged_decode_emissivity's.

broadband_emissivity(ε_vector, coefficients)

Collapse the five decoded narrowband emissivities (ASTER bands 10–14) to one broadband emissivity, the dot product with coefficients. NaN in any band propagates.

broadband_map(decoded_bands, coefficients)

Array form of broadband_emissivity and broadband_uncertainty: collapse a decoded (5, Nx, Ny) array (band index first, as in the HDF5 /Emissivity/* layout) to a broadband (Nx, Ny) array.

broadband_uncertainty(σ_vector, coefficients)

Propagate the five per-band standard deviations to a broadband uncertainty as the fully-correlated upper bound σ = Σ cᵢ σᵢ. ASTER TES band emissivities share one temperature retrieval and atmospheric correction, so their errors are strongly correlated.

fill_land_gaps(data, water, metadata)

Fill the clear-sky retrieval gaps over land from land alone. Water is blanked to NaN first so it cannot donate: a coastal land gap would otherwise inherit the open-ocean ε ≈ 0.98 against 0.94–0.96 over land. Water keeps its own retrievals, and its gaps stay NaN for the ungated inpainting that Field runs next.

place_tile!(field, tile_values, tile_longitude, tile_latitude, longitude, latitude)

Block-copy a decoded tile onto the regional grid (longitude, latitude) by mapping each tile-cell center to its native index. Indexing is by value, so the tile may store latitude in either direction; longitudes are folded into [-180, 180], so a tile lands on a grid labeled in any convention (this assumes the region spans less than 180° of longitude). NaN tile cells are skipped, leaving a valid value from an adjacent tile at a shared boundary in place, and cells outside the tile's footprint are untouched.

retrieve_data(metadata::ASTERGEDMetadatum)

Read the regional (Nx, Ny) broadband field for metadata.name and fill its land gaps from land. Water gaps are left NaN for the downstream inpainting. The decode and broadband collapse happen at download time, in the ArchGDAL extension.

albedo_read_window(metadatum)

Global native cell window (icols, jrows) (columns and ascending-frame rows) covering the metadatum's BoundingBox, mirroring construct_native_grid's restrict. Returns nothing for the global path (no region, or a longitude window that spans/wraps the ±180 seam), in which case the caller reads the whole global grid.

bluesky_blend(α_bs, α_ws, f)

Blend the black-sky (directional-hemispherical) albedo α_bs and the white-sky (bi-hemispherical) albedo α_ws into the blue-sky (actual) albedo, with f the diffuse fraction of the downwelling shortwave radiation:

\[α = (1 - f) α_bs + f α_ws\]

build_albedo_request(name, dates) -> Dict{String, Any}
build_albedo_request(
    name,
    dates,
    source_variables
) -> Dict{String, Any}

Construct the CDS request for the 1 km v2 albedo source_variables covering dates, which must share a (year, month) (CDS interprets year/month/nominal_day as a Cartesian product, and the valid nominal days differ by month). source_variables defaults to the full black-sky/white-sky pair, but download_ten_day_albedo! requests one variable at a time — see the note there.

download_ten_day_albedo!(
    retrieve,
    metadata::Metadata{<:CopernicusAlbedo};
    skip_existing,
    cleanup
) -> Any

Download the ten-day black-sky and white-sky broadband albedo for every date of metadata through retrieve(request, path) (supplied by a CDS backend extension), and repack each date's pair into a single compact local NetCDF. Requires ~/.cdsapirc credentials and acceptance of the Copernicus Global Land product license on the CDS portal.

Issues one request per source variable per calendar month, rather than one request for the pair: CDS answers a multi-variable request with a ZIP of per-variable members, and a backend is free to unwrap it to a single file — CopernicusClimateDataStore keeps only the first member — which silently drops half of the pair and fails the repack below.

repack_albedo_pair(blacksky, whitesky, destination_names, filepath, expected_size)

Repack the broadband albedo bands from a downloaded black-sky/white-sky product pair into one compact local NetCDF at filepath. blacksky and whitesky are (path, name) tuples locating each source file and its albedo variable; destination_names are the canonical local names (AL_DH_BB, AL_BH_BB). The packed integer representation, CF decoding attributes, and north→south latitude order of the source are preserved.

conversion_units(metadatum::Metadatum{<:ECCODarwin})

Set up conversion from the ECCODarwin output data to standard units

• salinity = SALTanom + 35
• biogeochemical tracer concentrations are in uL => umol/L in the output files from Darwin

metadata_filename(dataset, name, date, region)

Generate the filename for a given ECCO Darwin dataset and date.

The filename is constructed using the dataset variable name, and the iteration number is calculated from the date and epoch.

PrescribedAtmosphere(bounding_box, dates, dataset::ERA5PressureLevelsDataset; kw...)

Dataset-dispatched constructor: build an ERA5PrescribedAtmosphere over bounding_box at dates on dataset's native grid. Keyword arguments flow to ERA5PrescribedAtmosphere.

FieldTimeSeries(metadata::ERA5PressureMetadata, grid::PressureLevelGrid; kw...)

Load an ERA5 pressure-level variable raw onto an existing native PressureLevelGrid (e.g. a parent ERA5PrescribedAtmosphere's grid), placing each pressure level straight into its k-slot — the same per-level load the atmosphere applies to its own fields. This lets an extra ERA5 field share the parent's grid so it composes in AbstractOperations with the atmosphere's state.

The generic FieldTimeSeries(metadata, grid) would instead interpolate at load, freezing the data to one reference-height frame and desyncing it from the grid's time-varying geopotential heights.

download(
    metadata::Metadata{<:Union{ERA5HourlyLand, ERA5MonthlyLand}, D} where D;
    kw...
) -> String

Download ERA5-Land data for each date in metadata, returning the paths of the yearly files (repeated for dates that share a year).

download(
    metadatum::Metadatum{<:Union{ERA5HourlyLand, ERA5MonthlyLand}};
    skip_existing,
    kw...
) -> String

Return the path of the yearly ERA5-Land file containing metadatum, downloading it via download_era5_land unless skip_existing (default true) finds it already on disk.

batch_datetimes_for_cds(
    datetimes,
    dataset,
    num_vars;
    max_fields
) -> Vector{Vector{DateTime}}

Split datetimes into contiguous batches that each fit in one CDS request: each batch shares a (year, month) and contains at most max_dts_per_cds_request(dataset, num_vars; max_fields) datetimes. Batches are returned sorted by their first datetime so the caller can iterate in chronological order.

download_era5_land(metadatum; kw...)

Download the yearly ERA5-Land file containing metadatum from the Copernicus Climate Data Store and return its path.

Implemented in ext/NumericalEarthCopernicusClimateDataStoreExt.jl when CopernicusClimateDataStore is loaded; the fallback below fires only when the extension is not active.

foreach_nc(f, download_path, cleanup_dir) -> Any

If download_path is a ZIP archive (as CDS returns when mixing variable step types), extract all NetCDF files and call f(nc_path) on each. Otherwise call f directly on download_path.

group_by_calendar_month(datetimes) -> Dict

Group an iterable of DateTimes by (year, month). Returns a Dict whose keys are Tuple{Int, Int} (year, month) pairs and whose values are the datetimes that fall in that month. The 00:00 instant of a day belongs to that day (not the previous one).

max_dts_per_cds_request(
    dataset,
    num_vars;
    max_fields
) -> Any

Maximum number of datetimes that can share a single CDS request before the per-request cost limit is hit. Pressure-level datasets multiply by the number of selected pressure levels; single-level datasets count as one level. Falls back to 1 when a single datetime already exceeds the cap (which would force the caller's single-datetime download path).

mean_geopotential_heights(metadata::ERA5PressureMetadata)

Compute spatially and temporally averaged geopotential heights (m) for each pressure level in metadata.

Downloads the :geopotential field for every date in metadata, divides by g, averages over the horizontal domain and all dates, and returns one representative height per pressure level in bottom-to-top order (k=1 is highest pressure).

per_column_geopotential_discretization(metadata::ERA5PressureMetadata;
                                       clock = Clock{eltype(metadata)}(time = 0))

Build a PressureLevelVerticalDiscretization whose per-column heights follow the ERA5 geopotential in time. The geopotential Φ(λ,φ,p,t) is loaded over the whole date window of metadata (all snapshots in memory) and wrapped in a TimeSeriesInterpolation bound to clock; znode reads Φ interpolated to clock.time and divides by g. Sub-surface levels are clipped to the local surface geopotential so that columns stay monotonic. The single-level surface geopotential (orography × g, static in time) is downloaded from the matching ERA5*SingleLevel dataset and used as the clip source.

A static-z copy of the pressure-level dataset backs the Φ FieldTimeSeries, to break the recursive Field(ϕ) → native_grid → z_interfaces chain. The static z values are arbitrary placeholders; only their shape matters.

Pass the atmosphere's clock so that stepping the atmosphere advances the grid's heights; the default (a detached clock at t = 0) yields a first-date-static grid.

read_era5_yearly_series(paths, requested_times, name)

Read requested_times for variable name from paths, which is either a single file path (all requested times live in one file) or a Vector with one path per requested time (consecutive equal paths — e.g. every hour of the same year — are read together in a single, contiguous-when-possible slice, so each underlying file is opened at most once).

Returns (raw, λc, φc) where raw has shape (Nx, Ny, Nt) with latitude still in on-disk order (not yet reversed to -90→90).

split_era5_nc(src_path, nc_name_path_pairs, coordinate_vars)

Split a multi-variable NetCDF into individual per-variable files (single time step).

split_era5_nc_by_datetime(
    src_path,
    nc_varname_dt_path_triples,
    coordinate_vars,
    time_dimnames
)

Split a multi-timestep NetCDF into individual per-variable, per-timestep files, locating each requested DateTime in the file's own time coordinate. nc_varname_dt_path_triples is a vector of (nc_varname, datetime, dst_path). Because CDS expands year/month/day/time into a Cartesian product, the file may contain more datetimes than were asked for; matching by time value (not by position) extracts exactly the requested ones, and errors loudly if a requested datetime is absent.

split_era5_nc_multistep(
    src_path::AbstractString,
    nc_varname_tidx_path_triples,
    coordinate_vars,
    time_dimnames
)

Split a multi-timestep NetCDF into individual per-variable, per-timestep files. nc_varname_tidx_path_triples is a vector of (nc_varname, time_index, dst_path); src is a path or an already-open NCDatasets.NCDataset.

build_filename(dataset::ERA5LandDataset, name, dates::AbstractArray, region)

One yearly file covers every date within that year, so a request spanning multiple years resolves to multiple (repeated) filenames — one per date, naming that date's year file. set! for the resulting FieldTimeSeries groups consecutive same-year dates together and opens each yearly file once (see read_era5_yearly_series in ERA5_field_time_series.jl).

retrieve_data(metadata::ERA5Metadatum)

Retrieve ERA5 data from NetCDF file according to metadata. ERA5 is 2D surface data, so we return a 2D array with an added singleton z-dimension.

retrieve_data(metadata::ERA5PressureMetadatum)

Retrieve ERA5 pressure-level data from a NetCDF file. Returns a 3D array (lon, lat, level) with levels ordered bottom-to-top (highest pressure at k=1, lowest pressure at k=Nz).

retrieve_data(metadatum::Metadatum{<:ERA5YearlyFileDataset})

Read a 2D slice from the yearly ERA5 NetCDF file corresponding to the metadatum's date. Opens the yearly file, finds the time index matching the date, and extracts that timestep.

The yearly file contains all 8760-8784 hours for one year. This function indexes into the time dimension to extract just the requested hour.

set!(fts::ERA5NetCDFFTSMultipleYears, backend=fts.backend)

Load multiple timesteps from the yearly ERA5 file(s) into the FieldTimeSeries. Reads a time range from each yearly file in one operation (efficient!); a FieldTimeSeries spanning more than one calendar year reads from the corresponding number of yearly files (see read_era5_yearly_series).

This is called by Oceananigans when a FieldTimeSeries needs to load new data from disk.

AbstractGHSLDataset

Supertype for the Global Human Settlement Layer (GHSL) R2023A built-up rasters (European Commission JRC). All are distributed in the World Mollweide projection (ESRI:54009) and so require a reprojection to EPSG:4326 in the read path. Open access, no authentication.

Data source: https://human-settlement.emergency.copernicus.eu

built_surface_to_fraction(surface, cell_area)

Convert a GHS-BUILT-S built-up surface (m² of buildings within a native cell) to the plan-area index λᵖ = surface / cell_area, clamped to [0, 1]. cell_area is the native cell area in m² (100 × 100 at 100 m, 10 × 10 at 10 m). GHSL declares its no-data as a positive sentinel (65535 at 100 m, 255 at 10 m) that the warp has already written as NaN; a non-finite or negative surface maps to NaN.

julia> using NumericalEarth.DataWrangling.GHSL: built_surface_to_fraction

julia> built_surface_to_fraction(2500.0, 10_000.0)   # 25% built
0.25

julia> built_surface_to_fraction(10_000.0, 10_000.0)  # fully built, clamped
1.0

julia> isnan(built_surface_to_fraction(-1.0, 10_000.0))
true

ghsl_regional_raster(metadatum)

Geometry of the regional raster materialized for metadatum, taken from native_grid so that the file the reprojection writes lands on the cells the read path indexes. Returns the region to reproject onto — the native grid's extent, which reaches up to a cell beyond the requested window, and so is also the window whose tiles are needed — along with the cell counts Nx, Ny and the cell-center longitude and latitude.

julia> using NumericalEarth.DataWrangling: BoundingBox, Metadatum

julia> using NumericalEarth.DataWrangling.GHSL: GHSBuiltH, ghsl_regional_raster

julia> region = BoundingBox(longitude = (-0.11, -0.07), latitude = (51.505, 51.525));

julia> raster = ghsl_regional_raster(Metadatum(:building_height; dataset = GHSBuiltH(), region));

julia> (raster.Nx, raster.Ny)
(46, 25)

ghsl_tile_index(longitude, latitude)

Return the GHSL R2023A (row, column) tile index (1-based, from the NW corner) whose 1000 km Mollweide cell contains (longitude, latitude).

julia> using NumericalEarth.DataWrangling.GHSL: ghsl_tile_index

julia> ghsl_tile_index(2.35, 48.85)   # Paris
(4, 19)

ghsl_tile_url(dataset, row, column)

Return the JRC open-data HTTPS URL of the (row, column) GHSL tile .zip for dataset (each archive holds one Mollweide GeoTIFF of the built-up variable).

ghsl_tile_urls(dataset, region::BoundingBox)

.zip URLs of the GHSL tiles for dataset intersecting region, ready to be downloaded, unzipped, and mosaicked/reprojected by GDAL.

ghsl_tiles_in_bbox(region::BoundingBox)

Return the sorted (row, column) GHSL tiles whose 1000 km Mollweide cells intersect region, from the window's Mollweide (x, y) bounding box: y is monotone in latitude and |x| is largest at the equator-most latitude, so the extremes lie among the window's corners and its equator-most edge. Exact for regional windows; a complete (never-missing) superset for any window.

longitude_latitude_to_mollweide(longitude, latitude)

Forward World-Mollweide (ESRI:54009) projection: map (longitude, latitude) in degrees (central meridian 0°) to (x, y) in meters. The auxiliary angle θ solves 2θ + sin 2θ = π sin(latitude) by Newton iteration.

mask_building_height(value)

Map a GHS-BUILT-H ANBH value (meters) to a masked Float64, keeping every physical height (including the legitimate non-built value 0). GHSL declares its no-data as a positive sentinel (255) that the Mollweide→EPSG:4326 warp has already written as NaN; this pass carries that gap through and defensively maps any negative to NaN.

julia> using NumericalEarth.DataWrangling.GHSL: mask_building_height

julia> mask_building_height(12.5)
12.5

julia> mask_building_height(0)
0.0

julia> isnan(mask_building_height(NaN))   # the warp writes the no-data gap as NaN
true

resolution_meters(resolution::GHSBuiltSResolution)

The native pixel size in meters of a GHSBuiltSResolution.

globfp3d_parse_tile_bounds(name)

Parse a 3D-GloBFP tile filename gridID_lon1_lat1_lon2_lat2_…, whose coordinates are the SW and NE corners, into (; gid, west, south, east, north). Returns nothing if it does not match.

morphometry_latitude_bands(target_grid, region, Δ, maximum_raster_cells, padding)

Split the rows of target_grid covered by region into contiguous latitude bands whose native rasters (cell size Δ degrees) hold at most maximum_raster_cells cells each (a band never narrower than one row). Returns a vector of (; rows, region): the target-grid row range and the BoundingBox to rasterize for it — the full region longitude, and the rows' latitude interval widened by padding native cells but clamped to the region latitude, so every band raster is a sub-window of the single-pass raster and band rows reduce identically to a single pass.

reduce_morphometry(height, longitudes, latitudes, target_grid)

Reduce the fine building-height raster height (m), on the regular grid of cell-center longitudes/latitudes (degrees), onto target_grid. Returns the arrays that building_morphometry wraps in Fields; empty target cells are 0.

Fine cells are binned by the target cell faces, so a stretched target_grid works too.

retrieve_data(metadatum::JRA55Metadatum)

Read the 2D slice from the JRA55 NetCDF file corresponding to metadatum's single date. RepeatYearJRA55 resolves the index from the position within all_dates (the file holds exactly 2920 entries for 1990 and aligns 1:1 with all_dates); MultiYearJRA55 resolves the index against the file's own time axis (one Gregorian-calendar file per year).

read_cog_window(source, bbox)

Read the bbox longitude/latitude window from a single-band EPSG:4326 cloud-optimized GeoTIFF source, decode raw integers to physical units (mask nodata → NaN, then apply the band scale/offset), and return (longitude, latitude, data) with ascending, cell-center coordinates (latitude south-to-north, per CF convention).

Implemented in ext/NumericalEarthArchGDALExt/openlandmap.jl when ArchGDAL is loaded; the fallback below fires only when the extension is not active.

Statistic

Enum corresponding to the various prediction layers provided by the SoilGrids 2.0 dataset. Since SoilGrids was produced using quantile regression, multiple output statistics are available that capture uncertainty in estimates for each variable. The default is Mean which represents the most likely value. The 5th (Q5), 50th (Q50), and 95% (Q90) percentiles/quantiles represent the lower, middle, and upper bound of the predictions for each variable. These could be used in sensitivity analyses to check how sensitive model simulations are to uncertainty in soil properties. Note that the spread of the distribution should be expected to be larger in data-sparse regions such as remote islands, deserts, and the Arctic.

woa_z_interfaces_from_centers(depth_centers)

Compute cell interfaces (negative z, bottom-first) from WOA standard depth centers (positive, surface-first).

ESA_WORLDCOVER_CLASS_NAMES

The ESA WorldCover legend: verbose class name → LCCS class code, in ascending code order. ESAWorldCover tabulates the 11 classes.

The step between codes is not uniform (…90, 95, 100), so the codes are enumerated rather than derived from a stride. 0 is the no-data code and is not a valid class.

ESA_WORLDCOVER_VEGETATED_CLASSES

The classes summed into the :vegetation_fraction variable.

Which classes count as vegetated is a modeling choice rather than a property of the product. Bare/sparse vegetation (60) and moss/lichen (100) are excluded here because a land model driven by this fraction treats them as bare ground. A different partition needs no new download: every class also has its own :<class>_fraction variable, so any vegetated set is the sum of the corresponding fields.

aggregate_blockwise(pixels, factor, reduction)

Reduce the fine 2-D pixels raster onto a coarse grid by an integer factor, applying reduction (e.g. class_counts or majority_class) to each non-overlapping factor × factor block. Integer-factor aggregation keeps the coarse-cell boundaries aligned with fine-pixel boundaries — no reprojection of the categorical field.

size(pixels) must be divisible by factor in both dimensions.

aggregate_landcover(pixels, factor; vegetated_classes = ESA_WORLDCOVER_VEGETATED_CLASSES)

Reduce the fine pixels raster onto the coarse grid by an INTEGER factor, returning (; landcover_class, vegetation_fraction, class_fractions): the majority class code, the vegetated area fraction f_veg, and a NamedTuple of per-class area fractions keyed by the verbose class names. Each block is counted once, with every product derived from the same class_counts.

class_counts(pixels)

Return the number of pixels of each legend class in pixels, as a NamedTuple keyed by the verbose class names. Codes outside the legend — including the no-data code 0 — are counted by no class, so the sum of the counts is the number of valid pixels.

Every aggregated product is derived from these counts.

class_fraction_variable_name(class_name)

The per-class area-fraction variable name for class_name (e.g. :cropland → :cropland_fraction).

class_fractions(pixels)

Return a NamedTuple of per-class area fractions (each in [0, 1]) over pixels, keyed by the verbose class names. The fractions sum to 1 over valid pixels, and to 0 when every pixel is no-data.

majority_class(pixels)

Return the class code covering the most of pixels, ignoring no-data. Ties break toward the smaller code, and an all-no-data block returns 0.

This is the aggregation used for the categorical :landcover_class product — class codes are counted, never averaged.

vegetation_fraction(pixels; vegetated_classes = ESA_WORLDCOVER_VEGETATED_CLASSES)

Return the area fraction of pixels belonging to vegetated_classes — the subgrid f_veg. See ESA_WORLDCOVER_VEGETATED_CLASSES for the default set and why it is a modeling choice.

worldcover_window(longitude_bounds, latitude_bounds, factor)

Global native-pixel index bounds (i₁, i₂, j₁, j₂) of the read window covering a region at aggregation factor. The window snaps to the global aggregated-cell lattice and is padded by one aggregated cell on every side, guaranteeing it is a strict superset of the native grid that construct_native_grid builds (whose center-bracketing can extend up to one aggregated cell past the region edges). The read-back maps grid cells to file cells by an integer offset with no interpolation, so a narrower window would misregister the field by a whole cell.

atmosphere_ocean_salinity_flux(esm::EarthSystemModel)

Return the atmosphere-ocean salinity flux (g/kg m s⁻¹) at the atmosphere-ocean interface in a coupled esm.

atmosphere_ocean_temperature_flux(esm::EarthSystemModel)

Return the atmosphere-ocean temperature flux (K m s⁻¹) at the atmosphere-ocean interface in a coupled esm.

frazil_temperature_flux(esm::EarthSystemModel)

Return the two-dimensional frazil temperature flux (K m s⁻¹) in a coupled esm.

net_ocean_salinity_flux(esm::EarthSystemModel)

Return the net salinity flux (g/kg m s⁻¹) at the ocean's surface in a coupled esm.

net_ocean_temperature_flux(esm::EarthSystemModel)

Return the net temperature flux (K m s⁻¹) at the ocean's surface in a coupled esm.

sea_ice_ocean_salinity_flux(esm::EarthSystemModel)

Return the sea ice-ocean salinity flux (g/kg m s⁻¹) at the sea ice-ocean interface in a coupled esm.

sea_ice_ocean_temperature_flux(esm::EarthSystemModel)

Return the sea ice-ocean temperature flux (K m s⁻¹) at the sea ice-ocean interface in a coupled esm.

adopt_clock(component, clock)

Return component reconfigured so that its time is tracked with the same time type as the authoritative EarthSystemModel clock. EarthSystemModel construction calls this on every component so their clocks cannot drift apart over long runs — e.g. Float32 and Float64 clocks accumulating Δt differently across thousands of days.

Behavior depends on how the component tracks time:

• the generic method leaves component untouched, so a component with its own clock representation (e.g. SpeedyWeather or Veros, which track time internally) only extends this method if it needs coercion;
• a Simulation, whose clock type is fixed by its grid, errors on a mismatch since it cannot be coerced;
• prescribed components that carry an Oceananigans Clock extend this method through reclock, which coerces the clock to the model time type and warns when the type actually changes.

component_model(component)

Return the bare component model from a wrapper. ESM components are sometimes passed as a bare model (e.g. Breeze.AtmosphereModel, Breeze.RadiativeTransferModel) and sometimes as a Simulation wrapping that model (e.g. Simulation{<:Breeze.AtmosphereModel}). Component-interface methods that need the underlying model — to reach for .grid, .velocities, boundary conditions, etc. — call component_model(x) so they can share one implementation between the wrapped and unwrapped forms. The default unwraps a Simulation; the identity fallback covers bare models.

components(model::ESM)

Return a named tuple of the non-nothing components of an Earth System model.

default_earth_system_clock(
    atmosphere::Oceananigans.Simulations.Simulation
) -> Oceananigans.TimeSteppers.Clock{_A, _B, _C, Int64, Int64} where {_A, _B, _C}

Return the coupled model's default clock. A Simulation atmosphere's clock type is fixed by its grid and cannot be coerced, so the coupled clock adopts it (e.g. a Float32 atmosphere gets a Float32 coupled clock). Otherwise the clock defaults to Float64: prescribed atmospheres carry a coercible clock, so adopt_clock reconciles them to the authoritative model clock rather than the other way around.

materialize_earth_system_radiation!(atmosphere, radiation)

Return atmosphere with any radiation skeletons populated against the coupled-model radiation. Atmospheric components that carry an internal radiation handle (e.g. Breeze's AtmosphereModel.radiation) overload this to alias their internal flux divergence onto radiation.flux_divergence, so the atmosphere's tendency machinery reads from the same field the coupled radiation writes. Default: no-op, returning atmosphere unchanged.

Called from inside the EarthSystemModel constructor, after ComponentInterfaces is built (so the radiation has already been bound to the interfaces' surface temperature — see materialize_earth_system_surface_temperature).

materialize_earth_system_surface_temperature(
    radiation,
    interfaces
) -> Any

Return radiation with its surface temperature bound to the coupled interfaces' diagnostic surface (skin) temperature — the interface-formulation field the atmosphere actually sees, which coincides with a land's prognostic temperature only for bulk formulations — when radiation supports late binding and carries none of its own. Radiation components that read a surface temperature (e.g. Breeze's RadiativeTransferModel) overload this. Default: no-op, returning radiation unchanged.

materialize_sea_ice(sea_ice, ocean)

Return a fully-materialized sea_ice component that has access to the ocean's grid. Called from inside the EarthSystemModel constructor before ComponentInterfaces is built. Default: no-op, returning sea_ice unchanged. Sea ice types that require grid-allocated diagnostic fields (e.g. FreezingLimitedOceanTemperature) overload this.

abstract type AbstractInterfaceState{FT}

Interface state carried through the similarity-theory fixed-point solver (compute_interface_state). Concrete subtypes share the iterated quantities — fluxes (u★, θ★, q★), velocities (u, v), temperature (the skin temperature), and specific_humidity (qˢ) — and differ only in the surface property each interface needs: salinity for air–sea, the land hydrology / energy state for air–land.

AirIceInterfaceState{FT}

Air–sea-ice interface state. Sublimation is over fresh ice, so it carries no salinity (the Ice-phase saturation involves none, and the melting-point salinity the skin-temperature solve needs comes from the interior state). The humidity scalar is therefore zero.

AirLandInterfaceState{FT, H, E}

Air–land interface state. In place of salinity it carries the land's hydrology and energy surface state (e.g. (saturation = 𝒮,) and (temperature = Tᵢ,)), from which the surface humidity models derive what they need — the moisture availability β, the reservoir temperature, etc. β is not stored: it is evaporation_efficiency(efficiency, saturation), computed by the formulation.

AirLandRadiationState{FT}

Air-land interface radiation state at one cell: Stefan–Boltzmann constant σ, surface albedo α, emissivity ϵ, downwelling shortwave ℐꜜˢʷ, and downwelling longwave ℐꜜˡʷ. Returned by air_land_interface_radiation_state and consumed by the air-land flux kernel and apply_air_land_radiative_fluxes!.

AirSeaInterfaceState{FT}

Air–sea (ocean and sea-ice) interface state. Carries salinity, used by ImpureSaturationSpecificHumidity for the Raoult reduction of saturation.

AtmosphereSurfaceFluxes{F}

Atmosphere↔surface turbulent flux container, shared by the atmosphere–ocean, atmosphere–land and atmosphere–sea-ice interfaces: the five fluxes, and the three characteristic scales of the Monin–Obukhov solve.

ComponentExchanger(component, exchange_grid)

Hold a regridder, a buffer of state fields, and an optional correction used to bring data from a component (radiation, atmosphere, land, ocean, sea ice) onto a shared exchange_grid, where atmosphere–ocean and atmosphere–sea-ice fluxes are computed.

The optional correction is an in-place post-regrid hook applied to state after each interpolate_state! (e.g. AltitudeCorrection on the atmosphere). When correction === nothing, the per-step correct_state! sweep is a no-op for this component.

EdsonMomentumStabilityFunction{FT}

A struct representing the momentum stability function detailed by Edson et al. (2013). The formulation hinges on the definition of three different functions: one for stable atmospheric conditions $(ζ > 0)$, named $ψₛ$ and two for unstable conditions, named $ψᵤ₁$ and $ψᵤ₂$. These stability functions are obtained by regression to experimental data.

The stability parameter for stable atmospheric conditions is defined as

\[\begin{align*} dζ &= \min(ζ_{\max}, A⁺ ζ) \\ ψ⁺ &= - B⁺ ζ⁺ - C⁺ (ζ⁺ - D⁺) \exp(- dζ) - C⁺ D⁺ \end{align*}\]

While the stability parameter for unstable atmospheric conditions is calculated as a function of the two individual stability functions as follows

\[\begin{align*} f⁻₁ &= (1 - A⁻ζ)^{1/4} \\ ψ⁻₁ &= (B⁻ / 2) \log[(1 + f⁻₁ + f⁻₁² + f⁻₁³) / B⁻] - √B⁻ \mathrm{atan}(f⁻₁) - C⁻ \\ \\ f⁻₂ &= ∛(1 - D⁻ζ) \\ ψ⁻₂ &= (E⁻ / 2) \log[(1 + f⁻₂ + f⁻₂²) / E⁻]- √E⁻ \mathrm{atan}[(1 + 2f⁻₂) / √E⁻] + F⁻ \\ \\ f &= ζ² / (1 + ζ²) \\ ψ⁻ &= (1 - f) ψ⁻₁ + f ψ⁻₂ \end{align*}\]

The superscripts $+$ and $-$ indicate if the parameter applies to the stability function for stable or unstable atmospheric conditions, respectively.

EdsonScalarStabilityFunction{FT}

A struct representing the scalar stability function detailed by Edson et al. (2013). The formulation hinges on the definition of two different functions: one for stable atmospheric conditions $(ζ > 0)$, named $ψ⁺$ and one for unstable conditions, named $ψ⁻$.

These stability functions are obtained by regression to experimental data.

The stability parameter for stable atmospheric conditions is defined as

\[\begin{align*} dζ &= \min(ζ_{\max}, A⁺ζ) \\ ψ⁺ &= - (1 + B⁺ ζ)^{C⁺} - B⁺ (ζ - D⁺) \exp( - dζ) - E⁺ \end{align*}\]

While the stability parameter for unstable atmospheric conditions is calculated as a function of the two individual stability functions as follows

\[\begin{align*} f⁻₁ &= √(1 - A⁻ζ) \\ ψ⁻₁ &= B⁻ \log[(1 + f⁻₁) / B⁻] + C⁻ \\ \\ f⁻₂ &= ∛(1 - D⁻ζ) \\ ψ⁻₂ &= (E⁻ / 2) \log[(1 + f⁻₂ + f⁻₂²) / E⁻] - √E⁻ \mathrm{atan}[(1 + 2f⁻₂) / √E⁻] + F⁻ \\ \\ f &= ζ² / (1 + ζ²) \\ ψ⁻ &= (1 - f) ψ⁻₁ + f ψ⁻₂ \end{align*}\]

The superscripts $+$ and $-$ indicate if the parameter applies to the stability function for stable or unstable atmospheric conditions, respectively.

ImpureSaturationSpecificHumidity(phase [, water_mole_fraction=1])

Return the formulation for computing specific humidity at an interface.

InterfaceFluxScales{FT}

The solved similarity-theory characteristic scales at an interface: friction velocity u★, temperature flux scale θ★, and specific-humidity flux scale q★. Shared by every interface-state type.

InterfaceVelocities{FT}

The interface velocity (u, v) — the ocean surface current, or zero over land.

LogarithmicSimilarityProfile()

Represent the classic Monin–Obukhov similarity profile, which finds that

\[ϕ(z) = Π(z) ϕ_★ / ϰ\]

where $ϰ$ is the Von Karman constant, $ϕ_★$ is the characteristic scale for $ϕ$, and $Π$ is the "similarity profile",

\[Π(h) = \log(h / ℓ) - ψ(h / L) + ψ(ℓ / L)\]

which is a logarithmic profile adjusted by the stability function $ψ$ and dependent on the Monin–Obukhov length $L$ and the roughness length $ℓ$.

The exchange fluxes depend on the relative velocity between the atmosphere and the interface

ReynoldsScalingFunction(FT = Oceananigans.defaults.FloatType; A = 5.85e-5, b = 0.72)

Empirical fit of the scalar roughness length with roughness Reynolds number R_★ = u_★ ℓu / ν.

\[ ℓs = A / R_★ ^ b\]

See equation (28) by Edson et al. (2013).

SeaIceOceanInterface{J, F, T, S, P}

Container for sea ice-ocean interface data including fluxes, formulation, and interface state.

Fields

• fluxes::J: SeaIceOceanFluxes containing interfaceheat, frazilheat, salt, xmomentum, ymomentum
• flux_formulation::F: heat flux formulation (IceBathHeatFlux or ThreeEquationHeatFlux)
• temperature::T: interface temperature field (ocean surface view or computed field)
• salinity::S: interface salinity field (ocean surface view or computed field)

StateExchanger(grid, radiation, atmosphere, land, ocean, sea_ice;
               atmosphere_correction = nothing)

Container for one ComponentExchanger per component. The grid is the shared exchange grid onto which each component's state is regridded each time step. Per-component post-regrid corrections live on each ComponentExchanger and run as a sweep in phase 1.5 of the time step (see correct_state!).

TemperatureDependentAirViscosity([FT = Oceananigans.defaults.FloatType;
                                  ℂ₀ = 1.326e-5,
                                  ℂ₁ = ℂ₀ * 6.542e-3,
                                  ℂ₂ = ℂ₀ * 8.301e-6,
                                  ℂ₃ = - ℂ₀ * 4.84e-9])

Construct a TemperatureDependentAirViscosity object that calculates the kinematic viscosity of air as

\[ℂ₀ + ℂ₁ T + ℂ₂ T^2 + ℂ₃ T^3\]

WindDependentWaveFormulation(FT = Oceananigans.defaults.FloatType;
                             Umax = 19, ℂ₁ = 0.0017, ℂ₂ = -0.005)

A gravity wave parameter based on the wind speed ΔU with the formula ℂ₁ * max(ΔU, Umax) + ℂ₂ as shown in (Edson (2013)'s)[@cite Edson2013] equation (13) and surrounding text.

The exchange fluxes depend on the atmosphere velocity but not the interface velocity

Compute turbulent fluxes between an atmosphere and an interface state using similarity theory

Compute turbulent fluxes between an atmosphere and an interface state using similarity theory

buoyancy_scale(θ★, q★, ℂᵃᵗ, Tₛ, qₛ, g)

Return the characteristic buoyancy scale b★ associated with the characteristic temperature θ★, specific humidity scale q★, surface temperature Tₛ, surface specific humidity qₛ, atmosphere thermodynamic parameters ℂᵃᵗ, and gravitational acceleration g.

The buoyancy scale is defined in terms of the interface buoyancy flux,

\[u★ b★ ≡ w'b',\]

where u_★ is the friction velocity. Using the definition of buoyancy for clear air without condensation, we find that

\[b★ = (g / 𝒯ₛ) [θ★ (1 + δ qₛ) + δ 𝒯ₛ q★] ,\]

where $𝒯ₛ$ is the virtual temperature at the surface, and $δ = Rᵛ / Rᵈ - 1$, where $Rᵛ$ is the molar mass of water vapor and $Rᵈ$ is the molar mass of dry air.

Note that the Monin–Obukhov characteristic length scale is defined in terms of $b★$ and additionally the Von Karman constant $ϰ$,

\[L★ = u★² / ϰ b★ .\]

Calculate the air viscosity based on the temperature θ in Celsius.

compute_interface_heat_flux(flux::IceBathHeatFlux, ocean_state, ice_state, liquidus, ocean_properties, ℰ, u★)

Compute the heat flux at the sea ice-ocean interface using bulk formulation. Returns (Q, Tb, Sb) where:

• Q > 0 means heat flux from ocean to ice (ocean cooling)
• Tb, Sb are the interface temperature and salinity

compute_interface_heat_flux(flux::ThreeEquationHeatFlux, ocean_state, ice_state, liquidus, ocean_properties, ℰ, u★)

Compute the heat flux at the sea ice-ocean interface using three-equation formulation. Dispatches to the appropriate solve_interface_conditions based on whether the flux has internal conductive flux or not.

Returns (Q, Tb, Sb) where:

• Q > 0 means heat flux from ocean to ice (ocean cooling)
• Tb, Sb are the interface temperature and salinity

get_friction_velocity(u★, i, j, grid, τˣ, τʸ, ρᵒᶜ)

Return the friction velocity at grid point (i, j).

For a constant friction velocity (u★::Number), returns the value directly. For MomentumBasedFrictionVelocity, computes $u_* = \sqrt{|\tau| / \rho_o}$ from momentum stresses.

iterate_interface_state(flux_formulation, Ψₛⁿ⁻¹, Ψₐ, Ψᵢ, Qᵣ, ℙₛ, ℙₐ, ℙᵢ)

Return the nth iterate of the interface state Ψₛⁿ computed according to the flux_formulation, given the interface state at the previous iterate Ψₛⁿ⁻¹, as well as the atmosphere state Ψₐ, the interior state Ψᵢ, downwelling radiation Qᵣ, and the interface, atmosphere, and interior properties ℙₛ, ℙₐ, and ℙᵢ.

local_atmosphere_state(
    i,
    j,
    atmosphere_state,
    atmosphere_properties
) -> NamedTuple{(:z, :u, :v, :T, :p, :q, :h_bℓ), <:NTuple{7, Any}}

Atmospheric state at cell (i, j): the surface-layer velocities, temperature, pressure and specific humidity, with the elevation z at which they are taken and the boundary-layer height h_bℓ.

monin_obukhov_length(u★, b★, ϰ)

Return the Monin–Obukhov length scale

\[L★ = u★² / ϰ b★ ,\]

given the friction velocity u★, the buoyancy scale b★ (see buoyancy_scale), and the Von Karman constant ϰ. A neutral interface (b★ = 0) has an infinite L★, typed like u★ so that the profiles stay in the working precision.

sea_ice_ocean_interface(grid, sea_ice, ocean, flux_formulation)

Construct a SeaIceOceanInterface with the specified flux formulation.

For IceBathHeatFlux, the interface temperature and salinity point to the ocean surface values. For ThreeEquationHeatFlux, dedicated fields are created to store the computed interface values.

Arguments

• grid: the computational grid
• sea_ice: sea ice simulation
• ocean: ocean simulation
• flux_formulation: heat flux formulation (IceBathHeatFlux or ThreeEquationHeatFlux)

solve_interface_conditions(flux::ThreeEquationHeatFlux, Tᵒᶜ, Sᵒᶜ, ice_state, αₕ, αₛ, u★, ℰ, ρᵒᶜ, cᵒᶜ, liquidus)

Solve the three-equation system for interface temperature, salinity, and melt rate.

The three equations are:

1. Heat balance: $ρᵒᶜ cᵒᶜ αₕ u★ (Tᵒᶜ - T★) + κ (Tˢⁱ - T★) = ℰ q$
2. Salt balance: $ρᵒᶜ αₛ u★ (Sᵒᶜ - S★) = q (S★ - Sˢⁱ)$
3. Freezing point: $T★ = Tₘ(S★)$

where κ = k/(h ℰ) is the conductive heat transfer coefficient (zero for NoInternalFluxTEF).

Arguments

• ice_state: NamedTuple with fields S, h, hc, ℵ, T (internal temperature)

Returns (T★, S★, q) where q is the melt rate (positive for melting).

store_interface_fluxes!(
    interface_fluxes,
    interface_temperature,
    interface_specific_humidity,
    i,
    j,
    Ψₛ,
    Ψₐ,
    ℂᵃᵗ,
    ℒ,
    Tₛ,
    interface_properties
)

Write the turbulent fluxes and characteristic scales carried by the converged interface state Ψₛ into interface_fluxes at cell (i, j), the interface temperature Tₛ into interface_temperature, and the converged interface specific humidity into interface_specific_humidity. ℒ is the latent heat of the phase change at this interface — vaporization over water and land, sublimation over ice.

The sign convention is + for cooling of the interface and - for heating.

thermodynamic_constants(atmosphere)

Return the physical constants the atmosphere-state corrections need — the gravitational acceleration and the dry-air gas constant Rᵈ (the latter from the atmosphere's own thermodynamics) — so corrections like AltitudeCorrection don't hard-code or duplicate them.

validate_flux_formulation(
    fluxes::SimilarityTheoryFluxes,
    grid
)

Check that any Field-valued roughness length or zero-plane displacement of flux_formulation is laid out so the flux kernel can read it per cell on grid, and throw an ArgumentError naming the offending field otherwise.

clip_subsurface!(geopotential, surface_geopotential)

Clip each column of geopotential so that values below the local surface geopotential are replaced by the surface value. Required to keep columns monotonically increasing in z for the column bisection in _fractional_indices.

Works for either a Field (3-D geopotential at a single time) or a TimeSeriesInterpolation wrapping a FieldTimeSeries of geopotential.

HYPRES_REGRESSION

The Wösten et al. (1999) continuous pedotransfer regressions fitted to HYPRES, the database of HYdraulic PRoperties of European Soils (their Table 5). Each field predicts one parameter:

There is no regression for θʳ, and Wösten's regression for ηᴷ (R² = 12 %) is not included; HYPRESPedotransfer carries both as constants.

Coefficients follow the hypres_predictors order, one row per predictor family. A zero marks a predictor the paper's subset selection dropped.

WEYNANTS_REGRESSION

The Weynants et al. (2009) pedotransfer coefficients at the full precision of the spsh R package's ptf.cW, which their Table 6 rounds to four decimals, with the units of the Weihermüller et al. (2017) erratum. Each field predicts one parameter:

θʳ is absent because Weynants fitted it and found it "not significantly different from 0".

ComponentExchanger(land::SlabLand, grid)

Expose the generic atmosphere-facing SlabLand state: skin temperature T and surface saturation. Aerodynamic roughness lengths belong to the atmosphere-land flux closure (atmosphere_land_fluxes), not the land state.

abstract type AbstractFrontalAreaEstimator end

Supertype of the estimators of the frontal-area index λᶠ from the plan-area index λᵖ and mean height h.

A dataset that supplies only λᵖ leaves λᶠ to be assumed from a roughness-element shape — the dominant uncertainty of the drag partition. Where λᶠ is measured per cell it should be used directly instead.

abstract type AbstractHeightDistribution end

Supertype of the treatments of the roughness-element height distribution within a cell. UniformHeight takes every element to be as tall as the mean; VariableHeight derives the spread σʰ and the tallest element hᵐᵃˣ from the mean and corrects the uniform-height (ℓᵐ, d) for them.

HYPRESRegression(; porosity, inverse_air_entry_head, pore_size_uniformity,
                   matching_point_conductivity)

Coefficients of the continuous pedotransfer regressions fitted to HYPRES, the database of HYdraulic PRoperties of European Soils. One tuple per predicted parameter, each ordered by the shared predictor basis. porosity predicts ν directly; the other three predict ln αᵃᵉ, ln(𝓃 - 1) and ln K₀, the transforms that enforce αᵃᵉ > 0, 𝓃 > 1 and K₀ > 0.

Construct one to swap in a different calibration of the same functional form; HYPRES_REGRESSION holds the published fit.

abstract type PedotransferFunction

A pedotransfer function maps basic soil properties (texture, bulk density) to the van Genuchten–Mualem hydraulic parameters. Concrete subtypes implement

soil_hydraulic_parameters(ptf, sand, silt, clay, bulk_density, depth = 0)
    -> (; porosity, residual_liquid_fraction, inverse_air_entry_head,
          pore_size_uniformity, matching_point_conductivity,
          pore_connectivity_exponent)

with sand, silt, clay mass fractions (kg/kg), bulk_density in kg/m³, depth the soil layer's depth below the surface (m, positive down), and outputs in model units (αᵃᵉ in m⁻¹, K₀ in m s⁻¹), keyed like the keyword arguments of the hydrology closures.

WeynantsRegression(; porosity, inverse_air_entry_head, pore_size_uniformity,
                     matching_point_conductivity, pore_connectivity_exponent)

Coefficients of the Weynants et al. (2009) pedotransfer regressions, one tuple per predicted parameter, each ordered by the shared predictor basis (1, C, S, S², D, OC) in clay %, sand %, bulk density g/cm³ and organic carbon %. porosity and pore_connectivity_exponent are predicted directly; the other three predict ln αᵃᵉ, ln(𝓃 - 1) and ln K₀.

Construct one to swap in a different calibration of the same functional form; WEYNANTS_REGRESSION holds the published fit.

Interpolate the land state (freshwater fluxes) onto the exchange grid.

Scatter each prescribed freshwater component onto coastal ocean cells, conserving volume.

update_net_fluxes!(coupled_model, land::SlabLand)

Consume atmosphere-land turbulent fluxes and populate the precipitation, evaporation, vapor_flux, surface_energy_flux, and liquid_precipitation_flux accumulators declared by the land closures.

• surface_energy_flux ← 𝒬ᵀ + 𝒬ᵛ, positive upward (out of the slab).
• precipitation ← rainfall + condensation, positive into the slab.
• evaporation ← positive part of upward vapor flux.
• vapor_flux ← signed Jᵛ, positive upward (consumed by VariablySaturatedHydrology and WaterCoupledEnergy).
• liquid_precipitation_flux ← rainfall as Pˡ, positive downward (consumed by VariablySaturatedHydrology).

Radiative contributions are added on top in apply_air_land_radiative_fluxes!.

apply_height_distribution(
    _::UniformHeight,
    ℓᵐ,
    d,
    λᵖ,
    h
) -> Tuple{Any, Any}

Apply the height distribution to the uniform-height parameters (ℓᵐ, d) of a cell with plan-area index λᵖ and mean height h.

apply_height_distribution(
    _::UniformHeight,
    ℓᵐ,
    d,
    λᵖ,
    h,
    σʰ,
    hᵐᵃˣ
) -> Tuple{Any, Any}

Apply the height distribution to the uniform-height parameters (ℓᵐ, d) of a cell whose height spread σʰ and tallest element hᵐᵃˣ are measured rather than regressed from h.

build_diagnostic_accumulators(grid, energy, hydrology)

Allocate the closure-owned diagnostics NamedTuple (e.g. deep_liquid_flux, surface_runoff, water_storage_tendency). One closure can write a key the other reads (hydrology publishes water_storage_tendency; energy reads it for the variable-heat-capacity correction). Closures that need no diagnostics declare an empty tuple via diagnostic_variables(closure) = ().

build_flux_accumulators(grid, energy, hydrology)

Allocate the flux/forcing accumulator NamedTuple for a SlabLand. The coupler writes into these every step; closures only read.

build_flux_routing(target_grid, flux_time_series;
                   maximum_search_radius = 5, spread_radius = 1.2,
                   maximum_spread_cells = nothing)

Route a component stored as a per-area mass flux (kg m⁻² s⁻¹) on coastal cells — the JRA55-do convention — onto target_grid. Mouths are the cells positive in any record of flux_time_series, weighted by their source-cell area. Remaining keyword arguments go to build_river_routing.

diagnostic_variables(energy::AbstractEnergyBalance) -> Tuple{Vararg{Symbol}}

Names of closure-owned diagnostic fields published into SlabLand.diagnostics. Default: no diagnostics.

diagnostic_variables(hydrology::AbstractHydrology) -> Tuple{Vararg{Symbol}}

Names of closure-owned diagnostic fields published into SlabLand.diagnostics. Default: no diagnostics.

drag_partition_roughness_ratio(λᶠ, dʰ, Cᴰ, ϰ, β) -> Any

Roughness ratio ℓᵐ/h = (1 − dʰ)·exp{ −[0.5·β·(Cᴰ/ϰ²)·(1 − dʰ)·λᶠ]^(−1/2) }, obtained by partitioning the surface stress between the element drag and the ground and matching the log law over the canopy (Macdonald et al. (1998)). Non-monotonic in λᵖ: ℓᵐ rises then falls (isolated → wake-interference → skimming) with the frontal area λᶠ.

flux_variables(energy::AbstractEnergyBalance) -> Tuple{Vararg{Symbol}}

Names of flux/forcing accumulator fields this closure consumes from SlabLand.fluxes. The coupler writes into these every time step.

frontal_area_index(_::IsotropicFrontalArea, λᵖ, h) -> Any

Frontal-area index λᶠ of plan-area index λᵖ and mean height h under estimator.

height_spread(v::VariableHeight, h) -> Any

Height standard deviation σʰ (m) of a cell of mean element height h, from the affine regression σʰ = s₁·h + s₀, floored at zero below its one-story zero crossing.

height_spread_roughness(ℓᵐ, λᵖ, h, σʰ, a₁, b₁, c₁) -> Any

Rescale the uniform-height roughness length ℓᵐ (m) for the height spread σʰ by b₁·Y² + c₁·Y + a₁, with Y = λᵖ·σʰ/h (Kanda et al. (2013)). Reduces to a₁·ℓᵐ for a height-homogeneous canopy (σʰ → 0).

hydraulic_conductivity(
    i,
    j,
    grid,
    c::VanGenuchtenConductivity,
    𝒮,
    T
) -> Any

Darcy hydraulic conductivity (m s⁻¹) of closure c at saturation 𝒮 and temperature T (K) in cell (i, j).

initial_diagnostic(energy::AbstractEnergyBalance, name::Symbol, grid)

Build the initial diagnostic Field for diagnostic variable name. Defaults to a zeroed CenterField.

initial_diagnostic(hydrology::AbstractHydrology, name::Symbol, grid)

Build the initial diagnostic Field for diagnostic variable name. Defaults to a zeroed CenterField.

initial_flux(energy::AbstractEnergyBalance, name::Symbol, grid)

Build the initial flux Field for flux variable name. Defaults to a zeroed CenterField.

layer_depths(z_interfaces) -> Any

Depth below the surface of each layer's midpoint (m, positive down), deepest-first to match layer_weights.

using NumericalEarth.Lands: layer_depths

layer_depths([-1.0, -0.5, 0.0])

# output
2-element Vector{Float64}:
 0.75
 0.25

layer_weights(z_interfaces, slab_depth) -> Any

Per-layer thicknesses (m) inside the soil column [-slab_depth, 0], deepest-first to match the dataset vertical axis. z_interfaces are the layer faces increasing upward (e.g. [-1.0, -0.6, -0.3, 0.0]); layer k spans [z_interfaces[k], z_interfaces[k+1]], and layers outside the column get zero weight.

using NumericalEarth.Lands: layer_weights

layer_weights([-1.0, -0.6, -0.3, 0.0], 0.3)

# output
3-element Vector{Float64}:
 0.0
 0.0
 0.3

logexpm1(x) -> Any

log(eˣ - 1) for x > 0, evaluated as x + log(1 - e⁻ˣ) to keep eˣ from overflowing.

matched_retention_parameters(
    θ¹,
    θ²,
    θʳ,
    ν,
    heads
) -> Tuple{Any, Any}

The (αᵃᵉ, 𝓃) whose van Genuchten curve passes through water contents θ¹ and θ² at the two positive suction heads, given θʳ and ν. Eliminating αᵃᵉ leaves one equation in 𝓃,

log[(𝒮¹^(-1/𝓂) - 1) / (𝒮²^(-1/𝓂) - 1)] = 𝓃 log(heads[1] / heads[2]),   𝓂 = 1 - 1/𝓃,

whose residual increases monotonically in 𝓃; it is bisected over 1.01 ≤ 𝓃 ≤ 12, and αᵃᵉ follows from the first constraint. A root outside that bracket, or a pair of water contents that determines none, returns NaN.

maximum_element_height(v::VariableHeight, h, σʰ) -> Any

Tallest element hᵐᵃˣ (m) of a cell of mean height h and height spread σʰ, from the power-law regression hᵐᵃˣ = m₁·(σʰ)^m₂, floored at σʰ + h so the displacement parameter X stays within its [0, 1] range.

maximum_height_displacement(
    λᵖ,
    h,
    σʰ,
    hᵐᵃˣ,
    a₀,
    b₀,
    c₀
) -> Any

Displacement height d (m) referenced to the tallest roughness element hᵐᵃˣ, given the mean height h and the height standard deviation σʰ: d/hᵐᵃˣ = c₀·X² + (a₀·λᵖ^b₀ − c₀)·X, with X = (σʰ + h)/hᵐᵃˣ in [0, 1] (Kanda et al. (2013)).

normalize_property(::Type{FT}, prop)

Normalize a user-provided property value to a consistent in-memory representation: convert scalar numerics to type FT, keep fields as-is. Unsupported input types raise a MethodError at construction time.

outlet_indices_from_mask(outlet_mask, grid)

Return (outlet_i, outlet_j, outlet_λ, outlet_φ) for the true cells of a 2-D outlet_mask on grid.

packing_displacement_ratio(λᵖ, A, dᵐᵃˣ) -> Any

Displacement ratio d/h = 1 + A^(−λᵖ)·(λᵖ − 1) of an obstacle array packed to plan-area index λᵖ, clamped to [0, dᵐᵃˣ]. Empirical fit of Macdonald et al. (1998), satisfying d/h ≥ λᵖ and running from 0 at λᵖ = 0 to 1 at λᵖ = 1. Capping below 1 holds the displacement under roof level, and keeps ℓᵐ off the zero it would reach where the packed array closes into a smooth surface again.

property_value(prop, i, j[, k=1])

Evaluate a land property at (i, j, k). Defined for Number, AbstractArray, and AbstractField; unsupported types raise a MethodError at compile time. We avoid a catch-all throwing fallback here because GPUCompiler cannot lower the throw(ArgumentError("…$(typeof(p))…")) codepath (it needs runtime Symbol construction), which trips GPU kernel compilation even when the fallback isn't reached.

AbstractArray covers the GPU-adapted case: Oceananigans Fields are adapted to their underlying OffsetArray{T, 3, CuDeviceArray{T, 3}} on the device, which is not <: AbstractField even though it is the same memory.

saturation(hydrology::AbstractHydrology, land) -> AbstractField

Return the moisture-availability factor β ∈ [0, 1] the atmosphere's latent-heat BC reads. May be a Field, ZeroField, any AbstractField, or a Number.

spread_target_cells(wet, λₒ, φₒ, maximum_degrees, spread_radius, maximum_cells)

The ocean cells a mouth at (λₒ, φₒ) discharges into, in two steps: the mouth lands on the nearest wet cell within maximum_degrees, and the footprint is the wet cells within spread_radius degrees of that landing cell, nearest first, capped at maximum_cells. The footprint is centered on the landing cell rather than on the mouth, so a mouth relocated far offshore — the Ob and Yenisei move 2-3° — still gets a full footprint. spread_radius = nothing spreads over the whole maximum_degrees reach. Empty when no wet cell lies within maximum_degrees.

temperature_tendency(i, j, grid, energy, prognostic, fluxes, diagnostics, time)

Pointwise $∂Tˡᵃ/∂t$ for the energy closure energy, in kernel form. Dispatched on the closure; reads the land's prognostic, fluxes, and diagnostics NamedTuples. Each slab-temperature closure adds a method.

time_step!(energy, land, Δt[, time])

Advance the energy-balance state by Δt. The closure reads the land's prognostic fields (land.temperature, land.water_storage) and the surface energy flux (land.fluxes.surface_energy_flux) as needed. The default is a no-op (used by closures with no prognostic state, e.g. prescribed surface temperature).

update_diagnostics!(energy, land)

Refresh any cached diagnostics owned by the energy closure. Called by the container's update_state! at the end of each step. Default is no-op.

van_genuchten_m(𝓃) -> Any

The second van Genuchten shape parameter, 𝓂 = 1 - 1/𝓃, under the Mualem (1976) restriction.

van_genuchten_saturation(scaled_suction, 𝓃) -> Any

Effective saturation 𝒮 = [1 + scaled_suction^𝓃]^(-𝓂), with scaled_suction = -αᵃᵉ Π ≥ 0 the dimensionless suction; the inverse of VanGenuchtenRetention's pressure_head.

viscosity_correction(v::WaterViscosity, T) -> Any

Factor by which WaterViscosity scales hydraulic conductivity at temperature T (K), unity at v.reference_temperature.

using NumericalEarth
using NumericalEarth.Lands: viscosity_correction

v = WaterViscosity()

## warm soil conducts more than cold soil
round.([viscosity_correction(v, T) for T in (275.0, 293.0, 310.0)], digits = 3)

# output
3-element Vector{Float64}:
 0.603
 1.0
 1.453

set!(land::SlabLand; T=nothing, M=nothing)

Set the slab's prognostic skin temperature T and water storage M and refresh diagnostics in one call (so saturation is consistent with M afterward). Either keyword can be omitted to leave that field untouched. Each value is anything Oceananigans.set! accepts — a Number, Field, AbstractOperation, function (λ, φ, z), or array — so the initial state can stay abstract, e.g. set!(land; T = ERA5_T2m[1] - Γ * (z_land - z_era5), M = 75).

time_step!(land::SlabLand, Δt)

Advance the slab by Δt. Each closure runs its own time_step!, then update_state! refreshes diagnostics. Prognostic halos are filled at the end so atmosphere kernels reading the surface state see consistent values.

Closure-invocation order: hydrology → energy. Hydrology runs first so the energy step reads the same step's water_storage_tendency (dMˡᵃ/dt) and updated Mˡᵃ. The conservative WaterCoupledEnergy closure pairs the advective energy carried by each mass flux against cˡ(Tˡᵃ − Tᵣ) dMˡᵃ/dt; that cancellation — and hence energy conservation under water exchange at the slab temperature — is exact only when both use the mass flux hydrology actually applied this step.

The clock is ticked first; subsequent closures see land.clock.time = t + Δt. Time-dependent property providers should therefore evaluate at the post-step time (consistent with Atmospheres and Radiations).

update_state!(land::SlabLand)

Refresh closure-owned diagnostics. Order is hydrology → energy: hydrology produces saturation (𝒮) before the energy closure assembles any heat-capacity diagnostic from the freshly updated water storage.

Interpolated(source)

Wrap a parent state source for use as the condition value in an NormalFlowBoundaryCondition. source can be:

• a FieldTimeSeries (prescribed parent — interpolated in space + time), or
• an AbstractField (prognostic parent — interpolated in space at the parent's current state, which advances via time_step!).

During model construction Oceananigans' standard boundary-condition regularization tags this with the boundary's dimension / side / field location; afterward getbc evaluates the source at the appropriate boundary-face node.

NormalFlowBoundaryCondition(Interpolated(source); scheme = …) works with any scheme that consults getbc for its exterior value (e.g. PerturbationAdvection).

MultipleFluxes(flux_field, additional_fluxes)

A boundary-condition callable (intended to be wrapped in a discrete-form FluxBoundaryCondition) that combines two contributions at a field's top boundary:

1. flux_field: a 2D Field that an external flux solver (e.g. the OMIP coupled atmosphere/sea-ice solver) writes into each step.

2. additional_fluxes: any object compatible with getbc that returns a flux in the top cell

This lets the coupled flux solver and an additional flux share the same top boundary condition without one clobbering the other.

compute_absorption_coefficient!(radiation, time)

Refresh the blue-green absorption coefficient of radiation from its chlorophyll at time. Uniform chlorophyll carries its coefficient as a scalar and there is nothing to refresh.

default_or_override(default::Default, alternative_default=default.value) = alternative_default
default_or_override(override, alternative_default) = override

Either return default.value, an alternative_default, or an override.

The purpose of this function is to help define constructors with "configuration-dependent" defaults. For example, the default bottom drag should be 0 for a single column model, but 0.003 for a global model. We therefore need a way to specify both the "normal" default 0.003 as well as the "alternative default" 0, all while respecting user input and changing this to a new value if specified.

hydrostatic_ocean_simulation(grid;
                             clock = Clock(grid),
                             stop_time = default_stop_time(grid, clock),
                             Δt = estimate_maximum_Δt(grid),
                             closure = default_ocean_closure(),
                             tracers = (:T, :S),
                             free_surface = default_free_surface(grid),
                             reference_density = 1020,
                             rotation_rate = default_planet_rotation_rate,
                             gravitational_acceleration = default_gravitational_acceleration,
                             bottom_drag_coefficient = Default(0.003),
                             bottom_drag_background_velocity = 0,
                             implicit_bottom_drag = true,
                             forcing = NamedTuple(),
                             additional_surface_fluxes = NamedTuple(),
                             biogeochemistry = nothing,
                             timestepper = :SplitRungeKutta3,
                             coriolis = Default(HydrostaticSphericalCoriolis(; rotation_rate)),
                             momentum_advection = WENOVectorInvariant(time_discretization = AdaptiveVerticallyImplicitDiscretization(cfl=0.5)),
                             tracer_advection = WENO(order=7, time_discretization = AdaptiveVerticallyImplicitDiscretization(cfl=0.5)),
                             equation_of_state = TEOS10EquationOfState(; reference_density),
                             boundary_conditions::NamedTuple = NamedTuple(),
                             radiative_forcing = default_radiative_forcing(grid),
                             river_routing = nothing,
                             river_mouth_diffusivity = 0.1,
                             river_mouth_mixing_depth = 10,
                             warn = true,
                             verbose = false)

Construct and return a hydrostatic ocean simulation tailored to grid. Called by ocean_simulation(grid; model=:hydrostatic, ...).

This function assembles an Oceananigans's HydrostaticFreeSurfaceModel with physically consistent defaults for advection, closures, the equation of state, surface fluxes, Coriolis, barotropic pressure–gradient forcing, boundary conditions, and optional biogeochemistry. It then wraps the model into an Oceananigans's Simulation with the specified timestepping options.

Behavior and automatic configuration

Coriolis

• On spherical grids, an Oceananigans.Coriolis.HydrostaticSphericalCoriolis object is used by default.
• On rectilinear grids, Coriolis force is disabled unless explicitly provided.

Single-column grids (grid.Nx == 1 && grid.Ny == 1)

• Advection is turned off (momentum_advection = nothing, tracer_advection = nothing).
• Users may override bottom_drag_coefficient, but its default is 0.
• Immersed boundaries are ignored.

Bottom drag and immersed boundaries

For multi-column grids:

• Quadratic bottom drag is automatically applied to both u and v.
• Immersed-boundary bottom drag conditions are constructed for both velocity components.
• Barotropic potential forcings for u and v are also added automatically, and user forcing tuples (e.g. forcing = (u = ..., v = ...)) are appended if provided.

Radiative forcing

By default, radiative_forcing is TwoColorRadiation scheme.

River mouths

When river_routing is supplied, river_mouth_vertical_diffusivity is added to closure.

Tracers and closures

• tracers defaults to (:T, :S).
• If the closure requires turbulent kinetic energy (e.g. CATKEVerticalDiffusivity), the turbulent kinetic energy :e tracer is automatically added while its advection is disabled.

Boundary conditions

Default boundary conditions are constructed for u, v, T, and S, including surface fluxes and bottom drag. User-provided boundary conditions override the defaults on a per-field basis.

Keyword Arguments

• clock: Clock for the underlying model. Defaults to Clock(grid), a numeric clock starting at time = 0. Pass a DateTime-based clock to step the simulation in calendar time (e.g. when coupling).
• stop_time: Stop time for the simulation. Defaults to Inf for numeric clocks, or DateTime(9999, 12, 31, 23, 59, 59) for DateTime clocks. On Reactant architectures it defaults to nothing, since Reactant does not support stop_time.
• Δt: Timestep used by the Simulation. Defaults to the maximum stable timestep estimated from the grid.
• closure: A turbulence or mixing closure. Defaults to default_ocean_closure().
• tracers: Tuple of tracer names. Defaults to (:T, :S).
• free_surface: Free–surface solver. Defaults to default_free_surface(grid).
• reference_density: Reference seawater density used by the equation of state.
• rotation_rate: Planetary rotation rate used for Coriolis forcing.
• gravitational_acceleration: Gravitational acceleration, passed to buoyancy.
• bottom_drag_coefficient: Bottom drag coefficient. May be a Default wrapper.
• bottom_drag_background_velocity: unresolved velocity ub added in quadrature to the resolved speed in the quadratic drag, 𝛕 = - μ 𝐮 √(ub² + |𝐮|²), representing tides and other motions the grid does not carry.
• implicit_bottom_drag: whether the bottom and immersed quadratic drag are applied as affine fluxes with the drag coefficient carried in the vertical solver's diagonal. Default: true.
• forcing: Named tuple of additional forcing(s) for individual fields.
• additional_surface_fluxes: Named tuple of additional top boundary flux conditions (e.g. (; S=SurfaceFluxRestoring(...))) for any field (u, v, T, S).
• biogeochemistry: A biogeochemical model or nothing.
• timestepper: Time-stepping scheme; options are :SplitRungeKutta3 (default), or :QuasiAdamsBashforth2.
• coriolis: Coriolis object or Default(...) wrapper.
• momentum_advection: Momentum advection scheme. Defaults to WENOVectorInvariant(time_discretization = AdaptiveVerticallyImplicitDiscretization(cfl=0.5)).
• tracer_advection: Tracer advection scheme or named tuple of schemes. Defaults to WENO(order=7, time_discretization = AdaptiveVerticallyImplicitDiscretization(cfl=0.5)).
• equation_of_state: Equation of state object. Defaults to TEOS-10 (TEOS10EquationOfState).
• boundary_conditions: User-supplied boundary conditions; merged with defaults.
• radiative_forcing: Additional temperature forcing; merged into forcing.
• river_routing: NamedTuple of RiverRouting, typically land.river_routing. Defaults to nothing, which leaves closure untouched.
• river_mouth_diffusivity: vertical tracer diffusivity (m² s⁻¹) at the river mouths. Default: 0.1.
• river_mouth_mixing_depth: depth (m) over which it is applied. Default: 10.
• warn: If true, warnings are emitted for potentially unintended setups.
• verbose: If true, prints additional setup information.

merge_boundary_conditions(
    user::FieldBoundaryConditions,
    default::FieldBoundaryConditions
) -> FieldBoundaryConditions{_A, _B, _C, _D, _E, _F, _G, Nothing, Nothing} where {_A, _B, _C, _D, _E, _F, _G}

Merge user and default boundary conditions side-by-side: every side the user left unspecified (a DefaultBoundaryCondition) inherits the corresponding default side. This allows users to prescribe, for example, only the lateral boundary conditions of a field while retaining the default surface fluxes, bottom drag, and immersed boundary condition.

nonhydrostatic_ocean_simulation(grid;
                                clock = Clock(grid),
                                stop_time = default_stop_time(grid, clock),
                                Δt = 1,
                                closure = nothing,
                                tracers = (:T, :S),
                                coriolis = nothing,
                                reference_density = 1020,
                                gravitational_acceleration = default_gravitational_acceleration,
                                equation_of_state = TEOS10EquationOfState(; reference_density),
                                advection = WENO(order=9),
                                forcing = NamedTuple(),
                                boundary_conditions::NamedTuple = NamedTuple(),
                                verbose = false)

Construct and return a nonhydrostatic ocean simulation suitable for Large Eddy Simulation (LES). Called by ocean_simulation(grid; model=:nonhydrostatic, ...).

Uses Oceananigans' NonhydrostaticModel, which resolves the full 3D pressure field (no hydrostatic approximation). No free surface, barotropic forcing, or split-explicit timestepping is needed.

Keyword Arguments

• clock: Clock for the underlying model. Defaults to Clock(grid), a numeric clock starting at time = 0. Pass a DateTime-based clock to step the simulation in calendar time (e.g. when coupling).
• stop_time: Stop time for the simulation. Defaults to Inf for numeric clocks, or DateTime(9999, 12, 31, 23, 59, 59) for DateTime clocks. On Reactant architectures it defaults to nothing, since Reactant does not support stop_time.
• Δt: Timestep used by the Simulation. Defaults to 1 second.
• closure: Turbulence closure. Defaults to nothing (implicit LES via WENO advection scheme).
• tracers: Tuple of tracer names. Defaults to (:T, :S).
• coriolis: Coriolis parameter. Defaults to nothing.
• reference_density: Reference seawater density for the equation of state.
• gravitational_acceleration: Gravitational acceleration, passed to buoyancy.
• equation_of_state: Equation of state. Defaults to TEOS-10.
• advection: Advection scheme. Defaults to WENO(order=9).
• forcing: Named tuple of additional forcing(s).
• boundary_conditions: User-supplied boundary conditions; merged with defaults.
• verbose: If true, prints additional setup information.

allocate_interface_fluxes!(radiation, exchange_grid, surfaces)

Populate radiation.interface_fluxes with one InterfaceRadiationFlux per surface present in the model.

apply_air_land_radiative_fluxes!(coupled_model)

Add the radiative contribution to the land surface_energy_flux (positive upward) and write diagnostic radiative fluxes into coupled_model.radiation.interface_fluxes.land.

When coupled_model.radiation === nothing, this is a no-op.

apply_air_sea_ice_radiative_fluxes!(coupled_model)

Add the radiative contribution to the net sea-ice top heat flux and write the diagnostic radiative fluxes (upwelling LW, absorbed LW, transmitted SW) into coupled_model.radiation.interface_fluxes.sea_ice.

When coupled_model.radiation === nothing, this is a no-op. Also a no-op when sea-ice is not a Simulation{<:SeaIceModel}.

apply_air_sea_radiative_fluxes!(coupled_model)

Add the radiative contribution to the net ocean heat flux Jᵀ and write the diagnostic radiative fluxes (upwelling LW, absorbed LW, transmitted SW) into coupled_model.radiation.interface_fluxes.ocean.

When coupled_model.radiation === nothing, this is a no-op.

Interpolate the prescribed downwelling radiation onto the exchange grid.

set!(radiation::PrescribedRadiation; downwelling_shortwave=nothing, downwelling_longwave=nothing)

Set the prescribed downwelling shortwave and longwave radiative fluxes (W m⁻²), then refresh the interpolated state. Omitted keywords are left untouched. A Number sets a constant in space and time.

A thin proxy for the radiation interface that Breeze's AtmosphereModel expects. Holds an aliased flux_divergence field plus a back-reference to the source RadiativeTransferModel so that the atmosphere's update_radiation! call can delegate without the atmosphere knowing about the real RTM.

Defaults to the skeleton form (flux_divergence = nothing, radiative_transfer_model = nothing), which makes the atmosphere radiatively decoupled with zero extra allocation. materialize_earth_system_radiation! replaces the skeleton in-place when an EarthSystemModel is constructed with a non-nothing radiation.

atmosphere_model(grid;
                 surface_pressure = 101325,
                 potential_temperature = 285,
                 thermodynamic_constants = ThermodynamicConstants(eltype(grid)),
                 dynamics = CompressibleDynamics(; base_pressure = surface_pressure,
                                                 reference_potential_temperature = potential_temperature),
                 microphysics = SaturationAdjustment(equilibrium = WarmPhaseEquilibrium()),
                 momentum_advection = WENO(order=9),
                 scalar_advection = WENO(order=5),
                 boundary_conditions = NamedTuple(),
                 coriolis = nothing,
                 forcing = NamedTuple(),
                 closure = nothing,
                 clock = Clock{eltype(grid)}(time=0),
                 initialize = true)

Construct a Breeze AtmosphereModel with sensible defaults for coupled simulations. atmosphere_simulation wraps this in an Oceananigans Simulation (mirroring the role of ocean_simulation).

When initialize (the default) and dynamics isa CompressibleDynamics, the returned model is set to a resting, hydrostatically balanced state at the reference potential_temperature and surface_pressure, so its density is valid for anything that divides by it (e.g. the MOST surface-flux coupling). AnelasticDynamics needs no such step: its reference_state already prescribes a valid resting density, and zero prognostic perturbation is by construction the resting hydrostatically balanced state. Pass initialize = false when a caller derives the full state itself — the nested-child path does, via initialize_nested_child!, whose ERA5-derived reference state must not carry the resting initialization's default anchors.

Surface fluxes are handled by the EarthSystemModel coupling framework (via similarity theory), not by Breeze's own boundary conditions, so the bottom BCs on the prognostic momentum, energy, and moisture fields are pre-wired to 2D coupling fields that the coupler fills each step.

Radiation is wired in by the coupled-model constructor through materialize_earth_system_radiation!. The atmosphere is built with a skeleton CoupledRadiation proxy (no allocation, no tendency contribution), which is replaced inside EarthSystemModel with a materialized proxy that aliases coupled_model.radiation.flux_divergence. Passing a Breeze.RadiativeTransferModel directly as radiation here is rejected — use AtmosphereLandModel(atmosphere, land; radiation = rtm) instead.

To nest this child in a coarser parent atmosphere, use nested_atmosphere_model (Breeze extension), which derives the lateral BCs and Davies relaxation from the parent and wraps the result in a NestedModel.

atmosphere_simulation(
    grid;
    Δt,
    kw...
) -> SpeedyWeather.Simulation{V, M} where {V<:(Variables{Po, G, T, D, Pm, Pt, S, F} where {Po<:NamedTuple, G<:NamedTuple, T<:NamedTuple, D<:NamedTuple, Pm<:NamedTuple, Pt<:NamedTuple, S<:NamedTuple, F<:NamedTuple}), M<:PrimitiveWetModel}

Wrap atmosphere_model in an Oceananigans Simulation. Δt defaults to Inf because the coupled Simulation / NestedModel owns the time step in coupled use; pass a finite Δt to run! this Simulation directly. All other keyword arguments forward to atmosphere_model.

breeze_prognostic_state(
    constants::ThermodynamicConstants,
    pˢᵗ,
    T,
    qᵛ,
    qˡ,
    qⁱ,
    p
) -> NamedTuple{(:ρ, :θˡⁱ, :qᵗ), <:Tuple{Any, Oceananigans.Fields.Field{LX, LY, LZ, O, G, I, D, T, B, Oceananigans.Fields.FieldStatus{Float64}} where {LX, LY, LZ, O, G, I, D, T, B}, Any}}

Map a moist thermodynamic column state — temperature T, water-vapor specific humidity qᵛ, total liquid qˡ (cloud liquid + rain), total ice qⁱ (cloud ice + snow), and pressure p (all Fields on a common grid) — to the prognostic fields integrated by Breeze's CompressibleDynamics:

• density, ρ = p / (Rᵐ T), with Rᵐ = (1 − qᵗ) Rᵈ + qᵛ Rᵛ
• liquid-ice potential temperature, θˡⁱ = (T − (Lᵥ qˡ + Lₛ qⁱ) / cₚᵐ) / Π, with the Exner function Π = (p/pˢᵗ)^(Rᵐ/cₚᵐ) — Breeze's with_temperature, the exact inverse of the temperature the child reads it back with
• total-water specific humidity, qᵗ = qᵛ + qˡ + qⁱ

qˡ/qⁱ load the density through Breeze's mixture gas constant and set qᵗ (hence the dry density ρᵈ = ρ(1 − qᵗ)), so they must be the total liquid/ice mass fractions — every hydrometeor the parent carries, not just cloud condensate.

constants (Breeze.ThermodynamicConstants) supplies the gas constants, dry-air heat capacity, and latent heats; pˢᵗ is the potential-temperature reference pressure (the child dynamics' standard_pressure). Returns a NamedTuple (; ρ, θˡⁱ, qᵗ) of computed Fields. Both the ERA5 parent FieldTimeSeries population and the child initial condition in examples/era5_breeze.jl route through this so the conversion lives in one place.

nested_atmosphere_model(child_grid, parent_dataset; dates, kw...)

Build the parent PrescribedAtmosphere, nest a Breeze child in it, and initialize the child from parent_dataset at first(dates) — the returned model is ready to step. The parent spans child_grid's bounding box padded by parent_padding (default parent_dataset's default_horizontal_padding, margin for the lateral-BC interpolation stencils) at dates, on parent_dataset's native grid. Unless given, the default dynamics' surface_pressure anchor is the domain-mean dataset mean-sea-level pressure over the child at first(dates). When bottom_drag_coefficient is given, drag_surface_temperature defaults to the dataset's skin temperature at first(dates) regridded onto the child grid (a static snapshot, not the dataset's diurnal cycle). balancer controls the post-initialization adiabatic (DFI) balance: true (default) runs it, false skips it, and an AdiabaticBalancer(Δt=…) runs a custom (e.g. gentler) excursion. Remaining keyword arguments flow to nested_atmosphere_model(parent, child_grid; kw...).

nested_atmosphere_model(
    parent_atmosphere::PrescribedAtmosphere,
    child_grid;
    relaxation_rate,
    relaxation_width,
    relaxation_mask,
    sides,
    thermodynamic_constants,
    surface_pressure,
    reference_potential_temperature,
    terrain,
    terrain_blend_length,
    terrain_blend_width,
    terrain_smoothing_passes,
    bottom_drag_coefficient,
    drag_surface_temperature,
    parent_condensates,
    microphysics,
    momentum_advection,
    scalar_advection,
    coriolis,
    damping_rate,
    damping_depth,
    dynamics,
    boundary_conditions,
    forcing,
    kw...
) -> NestedModel{PrescribedAtmosphere{S, FT, G, T, U, Θ, Q, M, P, C, F, TP, TI}, AtmosphereModel{Dyn, Frm, Arc, Tst, Grd, Clk, Thm, Mom, Moi, Buy, Tmp, Sol, Vel, Trc, Adv, Cor, Frc, Mic, Cnd, Cls, Cfs, Rad, Prt}, NumericalEarthBreezeExt.StateExchanger{P1, Pr, C1, S1, Q1}} where {S, FT, G, T, U, Θ, Q, M, P, C, F, TP, TI, Dyn, Frm, Arc, Tst, Grd, Clk, Thm, Mom, Moi, Buy, Tmp, Sol, Vel, Trc, Adv, Cor, Frc, Mic, Cnd, Cls, Cfs, Rad, Prt, P1, Pr, C1, S1, Q1}

Build a Breeze child atmosphere over child_grid nested in parent_atmosphere, wrapped in a NestedModel. The child's prognostics (ρᵈ, ρu, ρv, ρθ, <moisture>) are precomputed from the parent's raw state on the parent grid as FieldTimeSeries (see child_prognostic_field_time_series); the child's lateral boundary conditions — and, when relaxation_rate (s⁻¹) is given, its interior Davies relaxation over relaxation_mask (default: a cosine ramp over the outermost relaxation_width cells) — interpolate those precomputed prognostics (via parent_boundary_conditions / parent_forcings).

Liquid/ice inputs to the combine default to the parent's hydrometeors — total liquid qᶜˡ + qʳ (cloud liquid + rain) and total ice qᶜⁱ + qˢ (cloud ice + snow) — but may be supplied from any source via parent_condensates, a NamedTuple with qᶜˡ/qʳ/qᶜⁱ/qˢ entries. Any missing or nothing entry — or the whole parent_condensates — is treated as absent (⇒ omitted; with all four absent, qᵗ = qᵛ).

Provides sensible, overridable physics defaults: microphysics (1-moment mixed-phase when CloudMicrophysics is loaded), momentum_advection = WENO(order=9), coriolis = SphericalCoriolis(), and a compressible split-explicit dynamics with an UpperSponge over the top damping_depth m at damping_rate; a matching ρw Rayleigh lid sponge (Relaxation toward zero) is added to forcing. Pass surface_pressure/reference_potential_temperature to anchor the default dynamics. Any boundary_conditions/forcing the caller passes are merged with the parent-derived ones (caller wins).

When bottom_drag_coefficient is given — a constant drag coefficient or a Breeze.PolynomialCoefficient — Breeze BulkDrag flux boundary conditions are applied at the bottom of the momentum densities ρu/ρv, computing the surface stress from the face-located near-surface velocity. drag_surface_temperature (a Field, function, or number) enters the drag's surface density and, for polynomial coefficients, its stability correction; compressible dynamics has no default surface temperature, so it must be supplied (the parent-dataset method defaults it to the dataset's skin temperature).

When terrain is given — an elevation Field, or a topography dataset (e.g. ETOPO2022()) that is regridded onto the child grid — the child grid's terrain-following coordinate is materialized in place before the model is built. The elevation is first smoothed with terrain_smoothing_passes applications of a binomial filter (smooth_topography!; 0 disables): point-sampled regridding otherwise leaves grid-scale orographic roughness that excites standing grid-scale noise in the near-surface flow on the terrain-following coordinate. If the parent knows its surface elevation (surface_elevation), the child elevation is then blended toward the parent's over an outer frame of physical width terrain_blend_length (meters; converted to a resolution-invariant cell count, or overridden directly with terrain_blend_width), so the terrain at the open boundaries matches the orography the parent state was produced with — and the blend slope stays fixed across resolutions rather than steepening.

atmosphere_simulation(spectral_grid::SpeedyWeather.SpectralGrid;
                      output_interval=nothing,
                      stop_time=nothing,
                      time_stepping=SpeedyWeather.Leapfrog(spectral_grid))

Return an atmosphere simulation using SpeedyWeather.PrimitiveWetModel on spectral_grid. Output is written when output_interval is provided. time_stepping controls SpeedyWeather's internal timestep (e.g. via its Δt_at_T32 field); the resulting Δt must be an integer divisor of the coupled EarthSystemModel timestep.

stop_time should match the stop_time later passed to the coupled Oceananigans.Simulation. It only synchronizes SpeedyWeather's internal clock (used for, e.g., progress.txt reporting); the actual run length is controlled by the coupled Simulation, not by SpeedyWeather's clock.

VerosOceanSimulation(setup, setup_name::Symbol)

Creates and initializes a preconfigured Veros ocean simulation using the specified setup module and setup name.

Arguments

• setup::AbstractString: The name of the Veros setup module to import (e.g., "global_4deg").
• setup_name::Symbol: The name of the setup class or function within the module to instantiate (e.g., :GlobalFourDegreeSetup).

install_veros()

Install the Veros ocean model Marine CLI using CondaPkg. Returns a NamedTuple containing package information if successful. Also patches Veros's signal handling to work with PythonCall.

patch_veros_signal_handling()

Patch Veros's signal handling to work with PythonCall. This prevents TypeError when Veros tries to set signal handlers from within Julia/PythonCall context. The issue is that Veros tries to set signal handlers, but PythonCall's wrapped Python objects aren't recognized as valid callable handlers by Python's signal.signal() function. We work around this by monkey-patching signal.signal() to accept any handler and convert invalid ones to SIG_DFL.

surface_grid(ocean::VerosOceanSimulation)

Constructs a LatitudeLongitudeGrid representing the surface grid of the given VerosOceanSimulation object. Notes: Veros always uses a LatitudeLongitudeGrid with 2 halos in both the latitude and longitude directions. Both latitude and longitude can be either stretched or uniform, depending on the setup, and while the meridional direction (latitude) is always Bounded, the zonal direction (longitude) can be either Periodic or Bounded.

Arguments

• ocean::VerosOceanSimulation: The ocean simulation object containing the grid state variables.

set!(ocean, v, x; path = :variable)

Set the v variable in the ocean model to the value of x. the path corresponds to the path inside the class where to locate the variable v to set. It can be either :variables or :settings.