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.

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.

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)

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 parameterises 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 or until the field is completely inpainted.

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, 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(metadata::Metadatum, grid::AbstractGrid; 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).

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.

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.

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.

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.

continue_downwards!(field, mask)

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

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

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

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, regions where mask[i, j, k] == 1 is inpainted and regions where mask[i, j, k] == 0 are preserved.

Arguments

• field: Field to be inpainted.

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

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

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.

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.

netrc_downloader(username, password, machine, dir)

Create a downloader that uses a netrc file to authenticate with the given machine. 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

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

retrieve_data(metadata)

Retrieve data from netcdf file according to metadata.

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.

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

conversion_units(metadatum::Metadatum{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}})

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.

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)

Build a PressureLevelVerticalDiscretization from the ERA5 instantaneous geopotential at the first date in metadata. The discretization stores raw Φ (m²/s²) and divides by g at read time. Sub-surface levels are clipped to the local surface geopotential so that columns are monotonic. The single-level surface geopotential is downloaded from the matching ERA5*SingleLevel dataset.

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

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

woa_z_interfaces_from_centers(depth_centers)

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

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_heat_flux(esm::EarthSystemModel)

Return the two-dimensional frazil heat flux (W m⁻²) 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.

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.

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 before ComponentInterfaces is built.

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 and atmosphere–land interfaces (both produce the same 8 quantities). Atmosphere–sea-ice uses a smaller container (AtmosphereSeaIceFluxes) because it does not emit the characteristic scales.

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. ElevationCorrection 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 and melt rate at the sea ice-ocean interface using bulk formulation. Returns (Q, q, Tᵦ, Sᵦ) where:

• Q > 0 means heat flux from ocean to ice (ocean cooling)
• q > 0 means melting (ice volume loss)
• Tᵦ, Sᵦ are the interface temperature and salinity

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

Compute the heat flux and melt rate 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, q, Tᵦ, Sᵦ) where:

• Q > 0 means heat flux from ocean to ice (ocean cooling)
• q > 0 means melting (ice volume loss)
• Tᵦ, Sᵦ 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 ℙᵢ.

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

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 ElevationCorrection don't hard-code or duplicate them.

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.

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.

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

update_net_fluxes!(coupled_model, land::SlabLand)

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

• net_energy_flux ← -(𝒬ᵀ + 𝒬ᵛ), positive into the slab.
• precipitation ← atmospheric rainfall flux + condensation (dew), both positive into the slab.
• evaporation ← net upward vapor flux, positive out of the slab.

Radiative contributions are added on top in apply_air_land_radiative_fluxes!.

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.

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.

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

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

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.

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.

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 net energy flux (land.fluxes.net_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.

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: energy → hydrology. Hydrology runs after energy so future closures that close the energy budget through phase change (snow melt, soil freeze/thaw) see the freshly updated temperature.

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.

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.

TwoColorRadiation(grid;
                  first_color_fraction,
                  first_absorption_coefficient,
                  second_absorption_coefficient)

Return TwoColorRadiation that computes the radiative flux divergence associated with a two-color radiation flux that decays according to Beer's law,

\[I(z) = ϵ₁ I₀ \exp(κ₁ z) + (1 - ϵ₁) I₀ \exp(κ₂ z)\]

where $I₀$ is the surface flux, $ϵ₁$ is the "first color" fraction of the total radiation, and $κ₁$ and $κ₂$ are the absorption coefficients for the two colors.

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;
                             Δ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),
                             forcing = NamedTuple(),
                             additional_surface_fluxes = NamedTuple(),
                             biogeochemistry = nothing,
                             timestepper = :SplitRungeKutta3,
                             coriolis = Default(HydrostaticSphericalCoriolis(; rotation_rate)),
                             momentum_advection = WENOVectorInvariant(),
                             tracer_advection = WENO(order=7),
                             equation_of_state = TEOS10EquationOfState(; reference_density),
                             boundary_conditions::NamedTuple = NamedTuple(),
                             radiative_forcing = default_radiative_forcing(grid),
                             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.

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

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

• Δ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.
• 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().
• tracer_advection: Tracer advection scheme or named tuple of schemes. Defaults to WENO(order=7).
• 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.
• warn: If true, warnings are emitted for potentially unintended setups.
• verbose: If true, prints additional setup information.

nonhydrostatic_ocean_simulation(grid;
                                Δ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

• Δ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 net energy flux 𝒬ˡᵃ 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.

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_simulation(grid;
                      surface_pressure = 101325,
                      potential_temperature = 285,
                      thermodynamic_constants = ThermodynamicConstants(eltype(grid)),
                      dynamics = AnelasticDynamics(ReferenceState(grid, thermodynamic_constants;
                                                                  surface_pressure, 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),
                      Δt = Inf)

Construct an Oceananigans Simulation wrapping a Breeze AtmosphereModel, with sensible defaults for coupled simulations. Mirrors the role of ocean_simulation.

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.

Returns the Simulation so callers can attach output writers, callbacks, or later wrap inside a coupled EarthSystemModel. The inner Δt defaults to Inf since the coupled Simulation (around an EarthSystemModel) owns the time step in coupled use; if you wrap this Simulation directly in a run!, pass a finite Δt.

atmosphere_simulation(spectral_grid::SpeedyWeather.SpectralGrid; output_interval=nothing)

Return an atmosphere simulation using SpeedyWeather.PrimitiveWetModel on spectral_grid. Output is written when output_interval is provided.

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.