Operators

Every physics process in AtmosTransport is implemented behind an abstract operator type with a No<Operator> no-op default. The runtime composes operators through a fixed Strang palindrome; users swap implementations by changing one field in the TOML config (or by implementing a new subtype and registering it).

The five operator families

The `apply!` contract

Each operator family's mutating entry point has a uniform shape, with the second positional argument carrying whatever time-resolved data that family needs:

Family	Signature
Advection	`apply!(state, fluxes, grid, scheme, dt; workspace, cfl_limit, diffusion_op, emissions_op, meteo)`
Diffusion	`apply!(state, meteo, grid, op::AbstractDiffusion, dt; workspace)`
Convection	`apply!(state, forcing::ConvectionForcing, grid, op::AbstractConvection, dt; workspace)`
Surface flux	`apply!(state, meteo, grid, op::AbstractSurfaceFluxOperator, dt; workspace = nothing)`
Chemistry	`apply!(state, meteo, grid, op::AbstractChemistryOperator, dt; workspace = nothing)`

Every method mutates state in place. Performance-sensitive implementations keep reusable scratch buffers in workspaces. Each family provides an identity operator (NoAdvection, NoDiffusion, NoConvection, NoSurfaceFlux, or NoChemistry) so the runtime can compose one typed model without conditionally calling absent physics.

Topology dispatch happens inside each apply! method, not on a separate axis: an LL state and a CS state route through different specialized kernels via Julia's multiple dispatch on the grid type.

Advection

AbstractAdvectionScheme is the root; the concrete schemes live in src/Operators/Advection/schemes.jl:

Subtype	Smooth-region reconstruction	Notes
`UpwindScheme`	First order	Donor-cell; cheap, very diffusive.
`SlopesScheme{L}`	Second order	Russell–Lerner face flux with the slope re-diagnosed from neighbouring means every sweep (van Leer). Not TM5's prognostic-moment slopes scheme. Limiter parameter `L`.
`PPMScheme{L,V}`	Second order	Limited piecewise-parabolic edges set a Russell–Lerner slope. Limiter parameter `L`. `V = FV3ScalarProfile` (CS only) integrates FV3's `kord = 8` parabola in the vertical sweep. Supported on LL and CS split-sweep; RG supports upwind only.
`LinRoodPPMScheme`	Piecewise parabolic	FV3 horizontal cross terms (CS only), paired by default with the FV3 `kord = 8` vertical profile in configured runs; explicit `vertical = "upwind"` retains first-order vertical transport. `ppm_order=5` or `7` selects the horizontal edge-value family; both share the same interior reconstruction.

Reconstruction order does not establish the temporal order or positivity of the complete transport update. See Advection schemes for the boundary treatment and Validation status for measured scheme coverage.

Limiter parameter L ranges over NoLimiter, MonotoneLimiter, PositivityLimiter and, for PPMScheme only, CW84Limiter (the complete Colella–Woodward PPM) — declared in the same file. PPMScheme() defaults to MonotoneLimiter(). The monotone and CW84 limiters are signed and constant-offset equivariant; only PositivityLimiter uses tracer zero as a bound.

TOML config (preferred form):

toml

[advection]
scheme = "slopes"     # "upwind" | "slopes" | "ppm" | "linrood"

# Cubed-sphere only: pick the Lin–Rood edge-value family (5 or 7).
# Only valid with scheme = "linrood"; setting ppm_order with
# scheme = "ppm" errors at config-parse time.
# scheme    = "linrood"
# ppm_order = 7

A NoAdvection identity operator is available for isolating other operators (e.g. convection-only timing, regression). Select with [advection] scheme = "none". Diffusion-only runs are supported as a single mass-conserving vertical solve. Surface fluxes require an advection palindrome and are therefore rejected with NoAdvection.

Diffusion

AbstractDiffusion is the root; concrete subtypes:

Subtype	Use
`NoDiffusion()`	Identity no-op; default when `[diffusion]` is absent or `kind = "none"`.
`ImplicitVerticalDiffusion{FT, KzF, SFC}`	Backward-Euler vertical diffusion driven by an `AbstractTimeVaryingField` Kz. `SFC` chooses midpoint source splitting or source-before-full-solve ordering.

Generic Kz diffusion uses a per-column backward-Euler tridiagonal solve for the new mixing ratio and updates tracer mass in flux form (each interface flux leaves one layer and enters the next), with a supported transpose in the CS adjoint. The exact TM5 Dkg path instead uses column-conservative bidiagonal factors directly on tracer storage, preserving weak transfers and signed mass budgets. Its reverse applies the factor transposes in reverse order. See Adjoint status for supported combinations.

TOML config ([diffusion] block):

toml

[diffusion]
kind  = "constant"
value = 1.0      # Kz [m²/s]; broadcast to all (i, j, k)

kind = "none" (or omitting the block entirely) selects NoDiffusion. Cubed-sphere binaries carrying PBL surface sections can use kind = "tm5_beljaars_viterbo_local_kz". On lat-lon and cubed-sphere grids, surface_flux_boundary = true adds configured surface mass before one full implicit vertical solve (S(dt) -> V(dt)) instead of using the separate midpoint split (V(dt/2) -> S(dt) -> V(dt/2)). The key name is historical: this changes operator ordering, not the lower row of the tridiagonal system. Reduced-Gaussian transport currently supports midpoint splitting only and rejects surface_flux_boundary = true while building the runtime recipe. Profile / derived / precomputed Kz fields exist in src/State/Fields/ — see State & basis for the full field-type list.

Convection

AbstractConvection is the root; concrete subtypes:

Subtype	Forcing carrier	Source
`NoConvection()`	—	Identity no-op; default.
`CMFMCConvection()`	`ConvectionForcing.{cmfmc, dtrain}`	GCHP-style upwind moist convection; mass flux + optional detrainment.
`CMFMCMatrixConvection()`	`ConvectionForcing.{cmfmc, dtrain}`	Implicit matrix solve assembled from GEOS cloud-mass flux and detrainment.
`TM5Convection{FT}()`	`ConvectionForcing.tm5_fields.{entu, detu, entd, detd}`	TM5 Tiedtke-1989 four-field entrainment / detrainment with an implicit column solve. Parametric on `FT`.

The convection subtypes consume a ConvectionForcing carrier (declared in src/MetDrivers/ConvectionForcing.jl) — different physics, identical plumbing. _refresh_forcing! populates model.convection_forcing at the first substep of each window by copying from the current met window; the operator does not call current_time itself.

TOML config ([convection] block):

toml

[convection]
kind = "cmfmc"     # or "cmfmc_matrix" / "tm5" / "none"

The runtime picks :cmfmc only against binaries whose header carries the :cmfmc payload section (and :dtrain if requested); similarly for :tm5 requiring :entu / :detu / :entd / :detd. Asking for a convection scheme the binary does not support is a load-time error, not a silent fallback.

Surface flux (sources)

AbstractSurfaceFluxOperator is the root; concrete subtypes:

Subtype	Use
`NoSurfaceFlux()`	Identity no-op; default.
`SurfaceFluxOperator{M}`	Applies a `PerTracerFluxMap` of `SurfaceFluxSource`s to the bottom-most layer (`k = Nz`).

SurfaceFluxOperator is materialized from [tracers.<name>.surface_flux] blocks and is the path for emissions inventories (EDGAR, GFED, GridFED, LMDz, …). See the worked CATRINE configs (config/runs/catrine_*.toml) for examples.

Inventories on the model grid

A static lat-lon inventory has flux densities $F_{s}$ [kg m⁻² s⁻¹] on source cells $s$ with areas $A_{s}$ : the file's cell_area or area variable, or, without one, the latitude-circle area $R^{2} Δ λ_{s} (\sin φ_{s + 1 / 2} - \sin φ_{s - 1 / 2})$ of the cell on the model sphere. The file's own total is $T = \sum_{s} F_{s} A_{s}$ ; for a per-cell inventory (EDGAR tonnes per cell) $F_{s} A_{s}$ is the stored amount itself. The source is the per-cell rate $r_{t}$ [kg s⁻¹] of model cell $t$ (surface_flux.regridding):

Conservative (default, every topology). The rate $m_{s} = F_{s} A_{s}$ is regridded as an extensive field with the intersection areas $W_{t s}$ of the source and model cells,
$r_{t} = \sum_{s} \frac{W_{t s}}{C_{s}} m_{s}, C_{s} = \sum_{t} W_{t s}, \sum_{t} r_{t} = \sum_{s} m_{s} = T .$
$C_{s}$ is the area of source cell $s$ that model cells cover (Regridding.source_coverage). The regridder's cells have great-circle edges, so a lat-lon cell's polygon is not its latitude-circle cell: regridding the density $F_{s}$ instead would conserve $\sum_{s} F_{s} A_{s}^{g c}$ , which differs from $T$ by up to $5 \times 10^{- 5}$ per 1° cell and $7 \times 10^{- 4}$ per 3.75°×1.875° cell (and by the file's radius convention, about $7 \times 10^{- 5}$ for the CAMS area). On reduced-Gaussian grids two rings with different longitude counts share a latitude as two different great-circle polylines, so a sliver between them belongs to no model cell and $C_{s}$ is smaller than the source cell (on O24, 6.3e-4 of the sphere is uncovered and a 1° cell loses up to 6.4 %; the density regrid gave a December 2021 CAMS total 1.2e-3 below the file's there). Dividing by $C_{s}$ gives each source cell's whole amount to the model cells it overlaps. On a target that does not cover the sphere (a regional or latitude-truncated mesh; the shipped runs are global) this also holds for a source cell that the target covers only in part: its whole amount goes into the covered part, whose density is raised by $A_{s} / C_{s}$ (+4.8 % in the boundary row of a 1° source on a lat-lon target ending at 74.75° N), and source cells that no model cell covers are dropped. Telling a sliver from a partial cover would need an area threshold, so this rule is kept and documented rather than special-cased.
Bilinear (lat-lon and reduced-Gaussian, only when set). The density is sampled at the model cell centre $x_{t}$ , ${\tilde{r}}_{t} = F (x_{t}) a_{t}$ , with the model cell area $a_{t}$ . The emission and uptake parts are renormalised separately,
$r_{t} = {\tilde{r}}_{t} \frac{P}{p} ({\tilde{r}}_{t} > 0), r_{t} = {\tilde{r}}_{t} \frac{N}{n} ({\tilde{r}}_{t} < 0), P = \sum_{s} max (m_{s}, 0), N = \sum_{s} min (m_{s}, 0),$
with $p$ and $n$ the same sums of ${\tilde{r}}_{t}$ . Both factors are positive, so no cell changes sign, and $\sum_{t} r_{t} = P + N = T$ . Renormalising the net total instead fails for fields with uptake: the factor $T / \sum_{t} {\tilde{r}}_{t}$ is arbitrary when the net is near zero and negative when the sampled net has the other sign, which turns every source into a sink (10 of 20 random fossil-plus-uptake fields on a 72×46 grid). When $p = 0 \neq P$ or $n = 0 \neq N$ no factor exists, and the build is an error; one uptake cell that falls between the sampled centres is enough. The factors keep the two totals, not the cell values, and can be far from 1: the uptake of the GridFED December 2021 file on a 5° LL72 grid has $N / n = 3.53$ , so every sampled uptake cell is multiplied by 3.53 (0.83 on 2.5°, 2.38 on 1°). Both factors are logged, and a factor more than 10 % from 1 is a warning. Point sources stay aliased: a 2° hotspot lands in whichever model cell's centre samples it.

Cubed-sphere grids always use the conservative path. Time-varying lat-lon series (cubed sphere only) are regridded slice by slice as densities with one regridder, so their slices conserve $\sum_{s} F_{s} A_{s}^{g c}$ ; their loader reads no file areas.

Time coverage of time-varying sources

A time-varying source evaluates its slices with clamped temporal schemes, so outside $[t_{1}, t_{n}]$ it holds its first or last slice. The loaders therefore require the run $[t_{s t a r t}, t_{e n d}]$ to lie inside the union of the slices' intervals, taken over all files of the source together. A slice with CF time bounds stands for its bounds, a months since stamp for its month (GriddedIO.known_slice_intervals). The stamps without bounds of all files are merged and sorted, and each holds for the regular step of that series (GriddedIO.stamp_intervals):

month starts step by one month, or by one year from a 1 January to the next: slice $k$ holds $[t_{k}, t_{k + 1})$ , the last slice one calendar month (one year after a yearly step);
other stamps step by their smallest spacing $Δ$ : slice $k$ holds $[t_{k}, t_{k + 1})$ when $t_{k + 1} - t_{k} = Δ$ (within 1 s), the last slice $Δ$ ;
a slice before a longer spacing (Jan 1 then Mar 1, or a missing day) holds for one step only, so the gap stays uncovered.

So daily files with one 00Z stamp each are a daily series, and only the very last stamp is extended. The period cannot be inferred, and the loaders stop with an error that asks for CF bounds, for a single stamp in all the files (a 1 January 00Z stamp may be a day, a month or a year) and for stamps about a month or a year apart that are not month starts (mid-month stamps: held from the stamp or centred on it?). A part of the run that no slice covers, including a gap between or inside files, is an ArgumentError naming the dates, unless the source sets hold_outside_coverage = true. The runner's run span ends 1 s after its last window for the rounding of its accumulated clock (RUN_SPAN_CLOCK_SLACK_SECONDS); the files need not cover that second.

Time-varying sources on a GPU

A TimeVaryingSurfaceFluxSource (cubed sphere; time_varying = true) stores its per-cell rates as slices $r_{c, k}$ at times $t_{k}$ . The series stays in host memory on every backend. On the device each source holds a small window buffer with the slices that the current met window $[t_{0}, t_{1}]$ reads, and the surface-flux kernels read the rates from it:

StepwiseFlux reads the slices in effect over the window, $k$ from $max {k : t_{k} \leq t_{0}}$ to $max {k : t_{k} \leq t_{1}}$ (clamped to the series).
LinearInterpFlux and ConservativeMeanFlux read from the left bracket slice at $t_{0}$ to the right bracket slice at $t_{1}$ .

$t_{1}$ is padded by $8 ε_{Float32} (t_{1} - t_{0}) + 4 ulp (t_{1})$ , because a Float32 step (for example 5400/7 s) can end slightly after the window. When a window becomes current, DrivenSimulation copies the slices of this range that the buffer does not hold yet, on the task that then launches the window's kernels. The buffer is a ring of slots (slice $k$ in slot $\mod 1 (k, capacity)$ ), so a slice that consecutive windows share is copied once, and its capacity is the largest window's slice count (two slices for stepwise, three for linear or conservative sources with daily or hourly series and hourly met windows), independent of the run length. The kernels see the same rates as with the whole series on the device, so results are bit-identical; a C90 2014–2024 daily series (0.78 GB per tracer in Float32) now costs a few 0.19 MB slices of device memory per tracer. The emission totals below sum the host series.

Emission totals

Sources hold per-cell model-storage rates $r_{c}$ (kg dry-air equivalent s⁻¹). Logs and budget diagnostics use global totals that the surface-flux module computes from the operator's own temporal reconstruction:

total_emission_rate(src, t): the global rate $R (t) = \sum_{c} r_{c} (t)$ at the instant t (seconds since the run start).
emission_segments(src, t, dt): the pieces $(R_{s}, τ_{s})$ in which one step $[t, t + d t]$ adds the source, with $\sum_{s} τ_{s} = d t$ .
integrated_emission(src, t0, t1, dt): the mass added over the steps $t_{n} = t_{0} + n d t$ , with $N = (t_{1} - t_{0}) / d t$ steps:

E = \sum_{n = 0}^{N - 1} \sum_{s} R_{s}^{(n)} τ_{s}^{(n)} .

A time-varying source stores slices with totals $R_{k} = \sum_{c} r_{c, k}$ . Every temporal scheme blends at most two slices, $R = w_{0} R_{i_{0}} + w_{1} R_{i_{1}}$ , so these totals equal the sum of the operator's per-cell increments up to summation order (test/core/test_surface_flux_integrals.jl checks this to 1e-14 relative). StepwiseFlux and ConservativeMeanFlux split each step at every slice time and evaluate each piece at its midpoint. This integrates the piecewise-constant or piecewise-linear reconstruction exactly, so their $E$ does not depend on $d t$ . LinearInterpFlux evaluates the rate at the step start, so its $E$ is a Riemann sum that depends on $d t$ .

The runner logs total_emission_rate(source, 0.0), the rate at the run start, for every source. Analysis.budget_closure (the series mode of scripts/diagnostics/tracer_budget_closure.jl), Analysis.carrier_drawdown and scripts/diagnostics/catrine_true_mass_balance.jl rebuild a run's sources and integrate them with these functions.

Chemistry

AbstractChemistryOperator is the root; concrete subtypes:

Subtype	Use
`NoChemistry()`	Identity no-op; default.
`ExponentialDecay`	Applies configured first-order loss rates to selected tracers.
`CompositeChemistry`	Applies a typed tuple of chemistry operators in sequence.

Constant-rate exponential decay uses exact half-steps before and after the transport/emissions window. Other chemistry, including exponential decay with a time-varying rate, runs after the transport and convection blocks. The same operator interface covers CellState and CubedSphereState; topology dispatch selects the storage-specific kernel launch.

Transport composition

The directional advection shell is a palindrome. Its center follows the diffusion operator's surface-coupling policy:

text

forward:   (X → Y → Z) × n   (fluxes / n; n = 1 unless CFL subcycling is needed)
center A:  V(dt/2) → S(dt) → V(dt/2)   (default symmetric Strang split)
center B:  S(dt) → V(dt)               (LL/CS source-before-full-solve policy)
           [without surface flux, the center is one V(dt); NoDiffusion is a no-op]
reverse:   (Z → Y → X) × n   (same n; the center runs once)
post:      apply!(convection)   (convection is outside the palindrome)
           apply!(chemistry)    (chemistry is outside the palindrome)

V is apply_vertical_diffusion_vmr! and S is apply_surface_flux!. Center A is time-symmetric and lets fresh bottom-layer storage diffuse in both half solves. Center B reproduces the supported GCHP-style placement of fresh surface storage before one full solve; it is an ordering policy, not a time-symmetric Strang center or a different tridiagonal boundary condition. Reduced-Gaussian transport implements Center A only.

Convection and chemistry sit outside the palindrome for the same reason: they are not commutative with advection at the per-substep level, and their natural cadence is the met window rather than the advection sub-step.

Adding a new operator

The same recipe applies in every family:

Subtype the abstract root: struct MyConvection <: AbstractConvection; … end.
Provide a No<Operator> peer if one doesn't already cover your slot — usually it does.
Implement apply!(state, …, op::MyConvection, dt; workspace) for whichever grid types you support. Multiple dispatch on the grid type handles topology specialization.
Wire selection from TOML through RuntimePhysicsSpecs.jl and its topology-dispatched materialize methods.
Test that the No<Operator> path is bit-exact to the explicit no-op — this is the contract that lets future code skip the slot for free.

What's next

Binary format — what the operator's input data looks like on disk, and how the runtime validates it.

Operators ​

The five operator families ​

The apply! contract ​

Advection ​

Diffusion ​

Convection ​

Surface flux (sources) ​

Inventories on the model grid ​

Time coverage of time-varying sources ​

Time-varying sources on a GPU ​

Emission totals ​

Chemistry ​

Transport composition ​

Adding a new operator ​

What's next ​

Operators

The five operator families

The `apply!` contract

Advection

Diffusion

Convection

Surface flux (sources)

Inventories on the model grid

Time coverage of time-varying sources

Time-varying sources on a GPU

Emission totals

Chemistry

Transport composition

Adding a new operator

What's next