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(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, use FieldTimeSeries(::MetadataSet).

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.

FieldTimeSeries(mset::MetadataSet, arch_or_grid=CPU(); kw...)

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

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.

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.

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

no_correction(i, j, grid, clock, model_fields, p)

Default correction function that returns zero. Used as the default value for the correction keyword arguments in os_papa_prescribed_flux_boundary_conditions.

ospapa_specific_humidity_fts(RHa, Ta, Pa, params)

Build a FieldTimeSeries of specific humidity (kg/kg) from OS Papa relative humidity (in %), air temperature (K), and pressure (Pa), using Thermodynamics.q_vap_from_RH with the Liquid() saturation curve.

vertical_interpolate(metadata::OSPapaMetadatum, z_src, data_src, z_dst)

Linearly interpolate a 1D profile from z_src centers onto z_dst levels. NaN values in data_src are skipped. Values outside the source range are extrapolated from the nearest valid value.

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.

components(model::ESM)

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

ComponentExchanger(component, exchange_grid)

Hold a regridder and a buffer of state fields 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.

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.

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)

Container for one ComponentExchanger per component. The grid is the shared exchange grid onto which each component's state is regridded each time step.

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

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.

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

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

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(),
                      kw...)

Construct a Breeze AtmosphereModel with sensible defaults for coupled simulations.

Surface fluxes are handled by the EarthSystemModel coupling framework (via similarity theory), not by Breeze's own boundary conditions.

Arguments

• grid: An Oceananigans grid for the atmosphere domain.

Keyword Arguments

• surface_pressure: Reference surface pressure in Pa. Default: 101325.
• potential_temperature: Reference potential temperature in K. Default: 285.
• thermodynamic_constants: Breeze ThermodynamicConstants. Default: ThermodynamicConstants(eltype(grid)).
• dynamics: Dynamics formulation. Default: AnelasticDynamics with the given reference state.
• microphysics: Microphysics scheme. Default: SaturationAdjustment(equilibrium = WarmPhaseEquilibrium()).
• momentum_advection: Advection scheme for momentum. Default: WENO(order=9).
• scalar_advection: Advection scheme for scalars. Default: WENO(order=5).
• boundary_conditions: Named tuple of boundary conditions. Default: NamedTuple().

All additional keyword arguments are forwarded to Breeze.AtmosphereModel.

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.