Operators API

For narrative coverage of the operator hierarchy, the apply! contract, and the Strang palindrome, see Operators. For per-scheme advection properties, Advection schemes.

`Operators` (top-level)

AtmosTransport.Operators.AbstractConvection Type

julia

AbstractConvection <: AbstractOperator

Root type for convective-transport operators. Concrete subtypes live in src/Operators/Convection/ (NoConvection, CMFMCConvection, TM5Convection).

AtmosTransport.Operators.AbstractDiffusion Type

julia

AbstractDiffusion <: AbstractOperator

Root type for vertical diffusion operators. Concrete subtypes live in src/Operators/Diffusion/operators.jl (NoDiffusion, ImplicitVerticalDiffusion, …).

AtmosTransport.Operators.AbstractOperator Type

julia

AbstractOperator

Root type for all physics operators in the transport model.

Advection

AtmosTransport.Operators.Advection.AbstractAdvectionScheme Type

julia

AbstractAdvectionScheme

Root abstract type for all advection operators in the mass-flux transport core.

Every concrete scheme belongs to one of three reconstruction families:

julia

AbstractAdvectionScheme
├── AbstractConstantScheme    (order 0: donor-cell / upwind)
├── AbstractLinearScheme      (order 1: van Leer slopes / MUSCL)
└── AbstractQuadraticScheme   (order 2: PPM / Prather moments)

This hierarchy enables orthogonal dispatch:

Reconstruction order selects the face-flux @inline function
Limiter type (carried as a type parameter) selects the slope/moment limiter
Grid topology is handled by the kernel shell (structured vs face-indexed)
CS execution style is handled separately (strang_split_cs! sweep shell vs Lin-Rood / FV3 horizontal update)
Backend (CPU/GPU) is handled by KernelAbstractions.jl

Implementing a new scheme

Subtype one of the three families
Implement _xface_tracer_flux, _yface_tracer_flux, _zface_tracer_flux (see reconstruction.jl)
The universal kernel shells in structured_kernels.jl will automatically dispatch to your face-flux functions at compile time

Example

julia

scheme = SlopesScheme(MonotoneLimiter())
strang_split!(state, fluxes, grid, scheme; workspace=ws)

AtmosTransport.Operators.Advection.AbstractConstantScheme Type

julia

AbstractConstantScheme <: AbstractAdvectionScheme

Piecewise-constant (order 0) reconstruction family.

The face value equals the donor cell mean — the cell upwind of the mass flux. This is the simplest conservative finite-volume scheme and the reference implementation for the generic kernel shells.

Concrete subtypes: UpwindScheme

AtmosTransport.Operators.Advection.AbstractLimiter Type

julia

AbstractLimiter

Policy object controlling slope and moment limiting in linear and quadratic advection schemes.

Limiters are carried as type parameters on scheme structs (e.g., SlopesScheme{MonotoneLimiter}), enabling compile-time specialization with zero runtime branches on GPU.

Available limiters:

NoLimiter: unlimited centered slopes (second-order, may oscillate)
MonotoneLimiter: van Leer minmod (monotone, TVD)
PositivityLimiter: ensures non-negative face values

See limiters.jl for the @inline implementations.

AtmosTransport.Operators.Advection.AbstractLinearScheme Type

julia

AbstractLinearScheme <: AbstractAdvectionScheme

Piecewise-linear (order 1) reconstruction family (van Leer 1977, MUSCL).

The subcell profile in each cell is $q (x) = \bar{q} + s_{x} (x - x_{c})$ where $s_{x}$ is a limited slope. The face flux is the Courant-fraction weighted integral of this profile over the swept volume (see _slopes_face_flux in reconstruction.jl).

Concrete subtypes: SlopesScheme

AtmosTransport.Operators.Advection.AbstractQuadraticScheme Type

julia

AbstractQuadraticScheme <: AbstractAdvectionScheme

Piecewise-quadratic (order 2) reconstruction family.

Includes PPM (Colella & Woodward 1984; Putman & Lin 2007) and Prather second-moment schemes. The subcell profile is a parabola constrained by the cell mean and (limited) edge values.

Concrete subtypes: PPMScheme (stub — kernels not yet implemented)

AtmosTransport.Operators.Advection.AdvectionWorkspace Type

julia

AdvectionWorkspace{FT, A, V1, A4}

Pre-allocated buffers for mass-flux Strang splitting. Eliminates all array allocations from the inner time-stepping loop.

The workspace provides TWO complete buffer pairs ((rm_A, m_A) and (rm_B, m_B)) used as ping-pong source/destination by strang_split!. Each directional sweep reads from one pair and writes to the other; the palindrome's six sweeps flip parity an even number of times, so the caller's home arrays receive the final result naturally.

Kernels always write to a DIFFERENT array than they read from — this is the double-buffer contract, and in-place updates would violate the stencil's read-before-write assumption and break mass conservation by ~10% per step.

Fields

rm_A::A, rm_B::A — 3D tracer-mass ping-pong pair (same size as rm)
m_A::A, m_B::A — 3D air-mass ping-pong pair (same size as m)
cluster_sizes::V1 — per-latitude clustering factors for reduced grids (Int32[Ny]; all ones for uniform grids; empty for face-indexed meshes)
face_left::V1, face_right::V1 — face connectivity for face-indexed meshes
rm_4d_A::A4, rm_4d_B::A4 — 4D tracer-mass ping-pong pair for the multi-tracer fused path ((Nx, Ny, Nz, Nt)). Both are allocated to size 0×0×0×0 when n_tracers == 0.
w_scratch::A — Thomas-factor scratch matching air_mass's shape: (Nx, Ny, Nz) for structured grids, (ncells, Nz) for face-indexed grids.
dz_scratch::A — layer-thickness input matching air_mass's shape. Caller is expected to fill this with current dz [m] (via hydrostatic from delp) before each diffusion apply!.

Constructors

julia

AdvectionWorkspace(m::AbstractArray{FT,3}; cluster_sizes_cpu=nothing, n_tracers=0)

Create workspace for a 3D structured grid, allocating both buffer pairs matching the size of m. If cluster_sizes_cpu is nothing, defaults to uniform (all ones).

julia

AdvectionWorkspace(m::AbstractArray{FT,2}; cluster_sizes_cpu=nothing, mesh=nothing)

Create workspace for a 2D face-indexed grid (cell × level layout).

AtmosTransport.Operators.Advection.AdvectionWorkspace Method

julia

AdvectionWorkspace(state::CellState; cluster_sizes_cpu=nothing, mesh=nothing)

Construct a workspace sized for state: infers n_tracers from ntracers(state) so the 4D ping-pong buffers match the packed tracer storage. This is the preferred form; the raw AdvectionWorkspace(m; n_tracers=…) is kept for low-level callers.

AtmosTransport.Operators.Advection.CSAdvectionWorkspace Type

julia

CSAdvectionWorkspace{FT, A, S, A4}

Pre-allocated cubed-sphere transport workspace.

rm_A, m_A are the halo-padded single-tracer advection ping-pong buffers shared across panels.
rm_4d_A is the packed-tracer panel buffer used by the production split-sweep path so CS follows the same packed tracers_raw paradigm as structured grids.
w_scratch, dz_scratch are panel-native column-operator workspaces with one structured (Nc, Nc, Nz) scratch array per panel.

AtmosTransport.Operators.Advection.LinRoodPPMScheme Type

julia

LinRoodPPMScheme{ORD} <: AbstractAdvectionScheme

Cubed-sphere Lin-Rood / FV3-style cross-term PPM advection with compile-time PPM order ORD.

This is distinct from PPMScheme: PPMScheme participates in the standard Strang split implemented by strang_split_cs!, while LinRoodPPMScheme selects the FV3-style horizontal Lin-Rood update (fv_tp_2d_cs!) paired with the existing vertical upwind sweep.

Supported orders currently match the implemented PPM edge-value families in ppm_subgrid_distributions.jl:

ORD = 5 — Huynh-constrained PPM
ORD = 7 — order-5 interior with special cubed-sphere face treatment

Examples

julia

LinRoodPPMScheme()    # default ORD=5
LinRoodPPMScheme(7)   # ORD=7 cubed-sphere boundary treatment

AtmosTransport.Operators.Advection.MonotoneLimiter Type

julia

MonotoneLimiter <: AbstractLimiter

Van Leer minmod limiter (van Leer 1977; Sweby 1984).

Limits the centered slope against the two one-sided differences scaled by 2, using the three-argument minmod function:

s = minmod (\frac{c_{i + 1} - c_{i - 1}}{2}, 2 (c_{i + 1} - c_{i}), 2 (c_{i} - c_{i - 1}))

This is TVD (total variation diminishing) and preserves monotonicity. The TM5 advectx__slopes / advecty__slopes routines use this limiter.

AtmosTransport.Operators.Advection.NoLimiter Type

julia

NoLimiter <: AbstractLimiter

No limiting applied. The slope is the full centered difference $s = (c_{i + 1} - c_{i - 1}) / 2$ . Second-order accurate but may produce new extrema (oscillations) near sharp gradients.

AtmosTransport.Operators.Advection.PPMScheme Type

julia

PPMScheme{L <: AbstractLimiter} <: AbstractQuadraticScheme

Piecewise Parabolic Method (Colella & Woodward 1984; Putman & Lin 2007).

Reconstructs a parabolic subcell profile constrained by the cell mean and limited edge values. Third-order accurate in smooth regions with appropriate limiting.

Status: structured-grid face-flux kernels are implemented and covered by kernel tests. PPMScheme is not yet part of the official real-data reference path, and face-connected support is still unimplemented.

Fields

limiter::L — parabolic profile limiting policy (default: MonotoneLimiter())

Example

julia

PPMScheme()                          # monotone-limited PPM
PPMScheme(NoLimiter())               # unlimited (high-order, may oscillate)

AtmosTransport.Operators.Advection.PositivityLimiter Type

julia

PositivityLimiter <: AbstractLimiter

Limits the slope to keep the reconstructed face values non-negative: $s = minmod (s, c_{i})$ , ensuring $c_{i} \pm s / 2 \geq 0$ .

Weaker than MonotoneLimiter but sufficient for species that must remain positive (e.g., tracer mixing ratios).

AtmosTransport.Operators.Advection.SlopesScheme Type

julia

SlopesScheme{L <: AbstractLimiter} <: AbstractLinearScheme

Van Leer / Russell–Lerner slopes advection (Russell & Lerner 1981).

Reconstructs a piecewise-linear subcell profile in each cell using a limited slope, then integrates the Courant-fraction swept volume to compute the face flux. Second-order accurate with MonotoneLimiter.

This is the method used by TM5 (advectx__slopes, advecty__slopes) for horizontal transport of atmospheric tracers.

The face tracer flux for positive mass flux $F > 0$ (left donor):

F_{q} = α (r_{m, L} + (1 - α) s_{x, L})

where $α = F / m_{L}$ is the Courant fraction, $r_{m, L}$ is the donor cell tracer mass, and $s_{x, L}$ is the limited first moment $s_{x} = m \cdot slope (c_{i - 1}, c_{i}, c_{i + 1})$ .

See _slopes_face_flux in reconstruction.jl for the full derivation and _limited_slope in limiters.jl for the limiter implementations.

Fields

limiter::L — slope/moment limiting policy (default: MonotoneLimiter())

Examples

julia

SlopesScheme()                       # monotone-limited (default, matches TM5)
SlopesScheme(NoLimiter())            # unlimited (2nd order, may oscillate)
SlopesScheme(PositivityLimiter())    # positivity-preserving

AtmosTransport.Operators.Advection.TracerView Type

julia

TracerView{FT, A}

Lightweight wrapper around a 4D array (Nx, Ny, Nz, Nt) that presents it as a 3D array (Nx, Ny, Nz) for a fixed tracer index t.

This enables the multi-tracer kernels to call the existing @inline face flux functions (which expect 3D rm arrays) without any code duplication. Julia's compiler inlines the getindex call and eliminates the wrapper entirely — zero overhead on both CPU and GPU.

Example

julia

rm_4d = zeros(Float64, 36, 18, 4, 50)  # Nx, Ny, Nz, Nt
rm_t = TracerView(rm_4d, Int32(3))       # view of tracer 3
rm_t[10, 5, 2]  # equivalent to rm_4d[10, 5, 2, 3]

AtmosTransport.Operators.Advection.UpwindScheme Type

julia

UpwindScheme <: AbstractConstantScheme

First-order donor-cell (Godunov) upwind scheme.

The tracer flux through a face is simply the mass flux times the mixing ratio of the upstream cell:

F_{q} = {\begin{cases} F \cdot c_{L} & if F \geq 0 \\ F \cdot c_{R} & if F < 0 \end{cases}

where $F$ is the mass flux [kg/s] and $c = r_{m} / m$ is the mixing ratio.

Properties: conservative, monotone, first-order accurate. Strongly diffusive but useful as a reference and for positivity-critical applications.

Example

julia

scheme = UpwindScheme()

AtmosTransport.Operators.Advection.apply_divergence_damping_cs! Method

julia

apply_divergence_damping_cs!(rm_panels, m_panels, mesh, ws, damp_coeff)

Conservative del-2 divergence damping on tracer panels. Mass-conserving flux-form Laplacian diffusion on mixing ratio (c = rm/m). Typical damp_coeff values: 0.02-0.05 for mild panel-boundary noise.

AtmosTransport.Operators.Advection.apply_linrood_update_adjoint! Method

julia

apply_linrood_update_adjoint!(lambda_rm, lambda_m,
                               lambda_fx_in, lambda_fx_out,
                               lambda_fy_in, lambda_fy_out,
                               lambda_rm_new, lambda_m_new,
                               am, bm, mesh)

Apply the discrete transpose of _linrood_update_kernel! for one panel: accumulate the adjoint of (rm_new, m_new) into the adjoint inputs (rm, m, fx_in, fx_out, fy_in, fy_out) for fixed velocities (am, bm).

All lambda_* adjoint accumulators are read-and-modified (atomically for face arrays) — callers are responsible for initialising them to zero before the call. The kernel only touches interior (i, j) indices 1..Nc; halo cells of lambda_rm/lambda_m and face cells outside 1..Nc+1 are left untouched.

AtmosTransport.Operators.Advection.apply_ppm_x_face_adjoint! Method

julia

apply_ppm_x_face_adjoint!(lambda_rm, lambda_m, lambda_fx_face, rm, m, am,
                            mesh, ::Val{ORD})

Discrete transpose of _ppm_x_face_kernel! (LinRood.jl:270) at ORD=5 or ORD=7. Folds _safe_mixing_ratio into the d6-AD chain and includes the donor-mass α = F / m_donor contribution. Atomic accumulation on shared cells.

ORD=7 applies the linear discontinuous-edge boundary correction at panel-edge faces (face_idx ∈ {1, Nc+1}) and recomputes the donor α-contribution against the corrected limited (q_L, q_R). Interior faces are bit-equal to ORD=5.

AtmosTransport.Operators.Advection.apply_ppm_x_face_from_q_adjoint! Method

julia

apply_ppm_x_face_from_q_adjoint!(lambda_q, lambda_fx_face, q, am, m,
                                   mesh, ::Val{ORD})

Discrete transpose of _ppm_x_face_from_q_kernel! (LinRood.jl:299) for one panel at ORD=5 or ORD=7 (LinRoodPPMScheme). The donor-mass denominator m_l/m_r in _ppm_face_value and the velocity am are treated as fixed parameters from the tape — the adjoint propagates only the q-stencil sensitivity. Atomic writes on lambda_q because multiple faces share each cell.

ORD=7 dispatches to a kernel that applies the linear discontinuous boundary correction (_apply_ord7_boundary_d6) at panel-edge faces (face_idx ∈ {1, Nc+1}) before the monotonicity step; interior faces are bit-equal to ORD=5.

AtmosTransport.Operators.Advection.apply_ppm_y_face_adjoint! Method

julia

apply_ppm_y_face_adjoint!(lambda_rm, lambda_m, lambda_fy_face, rm, m, bm,
                            mesh, ::Val{ORD})

Discrete transpose of _ppm_y_face_kernel! (LinRood.jl:241) at ORD=5 or ORD=7. See apply_ppm_x_face_adjoint! for the contract.

AtmosTransport.Operators.Advection.apply_ppm_y_face_from_q_adjoint! Method

julia

apply_ppm_y_face_from_q_adjoint!(lambda_q, lambda_fy_face, q, bm, m,
                                   mesh, ::Val{ORD})

Discrete transpose of _ppm_y_face_from_q_kernel! (LinRood.jl:325) for one panel at ORD=5 or ORD=7. See apply_ppm_x_face_from_q_adjoint! for the contract.

AtmosTransport.Operators.Advection.apply_pre_advect_x_adjoint! Method

julia

apply_pre_advect_x_adjoint!(lambda_rm, lambda_m, lambda_fx_face,
                              lambda_q_j, rm, m, am, fx_face, mesh)

Discrete transpose of _pre_advect_x_kernel! for one panel. See apply_pre_advect_y_adjoint! for the contract — same structure with X substituted for Y.

AtmosTransport.Operators.Advection.apply_pre_advect_y_adjoint! Method

julia

apply_pre_advect_y_adjoint!(lambda_rm, lambda_m, lambda_fy_face,
                              lambda_q_i, rm, m, bm, fy_face, mesh)

Discrete transpose of _pre_advect_y_kernel! for one panel: accumulate lambda_q_i into the adjoint accumulators of (rm, m, fy_face) for fixed velocity bm. The small-m_new zeroing exactly mirrors _safe_mixing_ratio (LinRood-style 100·eps threshold).

All lambda_* accumulators are read-modified; callers initialise them to zero before the call. Face writes use @atomic for shared neighbour-cell accumulation along Y.

AtmosTransport.Operators.Advection.copy_corners! Method

julia

copy_corners!(panels, mesh, dir)

Standalone corner fill for 6-panel fields. Rotates corner halo cells according to sweep direction dir (1 = X-sweep, 2 = Y-sweep).

AtmosTransport.Operators.Advection.diagnose_cm! Method

julia

diagnose_cm!(cm, am, bm, bt)

Backward-compatible wrapper for structured-grid continuity closure.

AtmosTransport.Operators.Advection.fill_panel_halos! Method

julia

fill_panel_halos!(panels::NTuple{6}, mesh::CubedSphereMesh; dir=0)

Fill halo regions of a 6-panel cubed-sphere field by copying interior data from neighboring panels with correct edge-to-edge orientation mapping.

Each panels[p] must be (Nc + 2Hp) × (Nc + 2Hp) × Nz with interior at [Hp+1:Hp+Nc, Hp+1:Hp+Nc, :].

If dir is 1 or 2, corner fill is performed for the given sweep direction (1=X, 2=Y) using the FV3 tp_core rotation formulas.

AtmosTransport.Operators.Advection.fillz_q! Method

julia

fillz_q!(q_panels, m_panels, mesh)

Post-advection positivity fixer for q-space transport. Fixes negative mixing ratios by borrowing mass from neighboring levels, then non-local column rescaling if needed. Port of GCHP's fillz (fv_fill.F90:51-156).

AtmosTransport.Operators.Advection.fv_tp_2d_cs! Method

julia

fv_tp_2d_cs!(rm_panels, m_panels, am_panels, bm_panels,
              mesh, ::Val{ORD}, ws, ws_lr; damp_coeff=0.0)

Lin-Rood horizontal advection for cubed-sphere grids. Averages X-first and Y-first PPM orderings (FV3 fv_tp_2d algorithm).

AtmosTransport.Operators.Advection.fv_tp_2d_cs_q! Method

julia

fv_tp_2d_cs_q!(q_panels, m_panels, am_panels, bm_panels,
                 mesh, ::Val{ORD}, ws, ws_lr; damp_coeff=0.0)

Q-space Lin-Rood horizontal advection. Evolves q (mixing ratio) and m (air mass, same role as pressure thickness) in-place. m_panels is both read (for CFL fraction in PPM face values) and written (mass divergence update).

AtmosTransport.Operators.Advection.reconstruction_order Method

julia

reconstruction_order(scheme) → Int

Return the polynomial order of the subcell reconstruction:

0 for constant (upwind)
1 for linear (slopes)
2 for quadratic (PPM)

Useful for diagnostics and for selecting stencil widths in multi-tracer kernel fusion.

AtmosTransport.Operators.Advection.required_halo_width Method

julia

required_halo_width(scheme) -> Int

Return the minimum cubed-sphere halo width needed by scheme's horizontal stencil. This is a capability query, not a reconstruction-order query: several schemes can share the same polynomial family while using different CS execution paths.

AtmosTransport.Operators.Advection.strang_split! Method

julia

strang_split!(state, fluxes, grid, scheme; workspace)

Perform one full Strang-split advection step on a structured mesh.

Splitting sequence

julia

  X → Y → Z → Z → Y → X
  ─────────────────────────
  half    half   full  half    half

All Nt = ntracers(state) tracers are advanced together in a single multi-tracer kernel launch per direction. The mass update is computed once per cell; tracer fluxes are evaluated per-tracer inside the kernel. The per-tracer Julia loop has been eliminated (plan 14 Commit 4) in favour of strang_split_mt! on the packed state.tracers_raw buffer.

Arguments

state::CellState — contains air_mass and tracers_raw
fluxes::StructuredFaceFluxState — mass fluxes (am, bm, cm)
grid::AtmosGrid{<:LatLonMesh} — structured lat-lon grid
scheme — advection scheme (AbstractAdvectionScheme)
workspace::AdvectionWorkspace — pre-allocated double buffers; use AdvectionWorkspace(state) so the 4D ping-pong buffers are sized for ntracers(state).

AtmosTransport.Operators.Advection.strang_split_cs! Method

julia

strang_split_cs!(panels_rm, panels_m, panels_am, panels_bm, panels_cm,
                 mesh, scheme, workspace; flux_scale=1, cfl_limit=0.95)

Perform one Strang-split advection step on a 6-panel cubed-sphere field with automatic CFL-based subcycling per direction.

Splitting sequence

julia

X sweep (n_x subcycles)
→ fill_panel_halos!(dir=1)     ← exchange halos between panels (X direction)
→ Y sweep (n_y subcycles)
→ fill_panel_halos!(dir=2)     ← exchange halos between panels (Y direction)
→ Z sweep (n_z subcycles)      ← first Z half-step
→ Z sweep (n_z subcycles)      ← second Z half-step (palindrome)
→ fill_panel_halos!(dir=2)
→ Y sweep (n_y subcycles)
→ fill_panel_halos!(dir=1)
→ X sweep (n_x subcycles)

This palindromic sequence (X-Y-Z-Z-Y-X) gives second-order temporal accuracy via Strang (1968) symmetry. The halo exchanges must happen BETWEEN successive horizontal sweeps because the panel-edge reconstruction stencil reads from adjacent panels.

Subcycling

Each direction D ∈ {X, Y, Z} has its own subcycle count n_D determined by an evolving-mass CFL pilot: the pilot applies n_D passes of flux_scale/n_D, checking that no cell mass goes negative or that |outgoing_flux| < cfl_limit × cell_mass at each pass. If the pilot fails, n_D is incremented until it passes (or hits max_n_sub and errors).

Panel array layout

Each panel's rm and m arrays are (Nc+2Hp, Nc+2Hp, Nz) with Hp-wide halos. Interior cells are at indices [Hp+1:Hp+Nc, Hp+1:Hp+Nc, :]. The sweep kernels only update interior cells; halo regions are filled by fill_panel_halos! from adjacent panels.

Arguments

panels_rm, panels_m: NTuple{6} of 3D arrays (Nc+2Hp, Nc+2Hp, Nz) — tracer mass and air mass. Modified in-place.
panels_am, panels_bm, panels_cm: NTuple{6} of flux arrays. am[Nc+2Hp+1, Nc+2Hp, Nz], bm[Nc+2Hp, Nc+2Hp+1, Nz], cm[Nc+2Hp, Nc+2Hp, Nz+1]. Read-only.
mesh: CubedSphereMesh with Nc, Hp, and panel connectivity.
scheme: advection scheme — UpwindScheme() uses gamma-clamped upwind (positivity-safe). SlopesScheme() and PPMScheme() use the generic KA kernels with _xface_tracer_flux dispatch. Higher-order schemes require mesh.Hp ≥ 2 (Slopes) or mesh.Hp ≥ 3 (PPM).
workspace: pre-allocated CSAdvectionWorkspace buffers.
flux_scale: overall scaling applied to all fluxes (default 1.0).
cfl_limit: maximum CFL per subcycle pass (default 0.95).

AtmosTransport.Operators.Advection.strang_split_cs_mt! Method

julia

strang_split_cs_mt!(panels_rm_4d, panels_m, panels_am, panels_bm, panels_cm,
                    mesh, scheme, workspace; ...)

Packed-tracer cubed-sphere split-sweep transport. This is the production CS path for CSSplitSweepStyle schemes: air mass is advanced once per sweep and all tracers in each panel's fourth dimension are updated inside the same panel kernel. The sequence and CFL contract match strang_split_cs!.

AtmosTransport.Operators.Advection.strang_split_linrood_ppm! Method

julia

strang_split_linrood_ppm!(rm_panels, m_panels, am_panels, bm_panels, cm_panels,
                           mesh, ::Val{ORD}, ws, ws_lr; cfl_limit=0.95, damp_coeff=0.0)

Full 3D advection: Horizontal(LR) → Z → Z → Horizontal(LR).

AtmosTransport.Operators.Advection.strang_split_mt! Method

julia

strang_split_mt!(rm_4d, m, am, bm, cm, scheme, ws)

Multi-tracer Strang-split advection on a packed 4D tracer array.

This is the performance-optimized path: all Nt = size(rm_4d, 4) tracers are processed in a SINGLE kernel launch per sweep direction (6 total), rather than 6 × Nt launches in the per-tracer path.

The mass update (m_new = m + flux_in - flux_out) is computed ONCE per cell per sweep, shared across all tracers.

Arguments

rm_4d::AbstractArray{FT,4} — tracer mass (Nx, Ny, Nz, Nt), mutated
m::AbstractArray{FT,3} — air mass (Nx, Ny, Nz), mutated
am, bm, cm — mass fluxes (x, y, z directions)
scheme::AbstractAdvectionScheme — advection scheme
ws::AdvectionWorkspace — workspace with 4D buffer allocated

Convection

AtmosTransport.Operators.Convection.CMFMCConvection Type

julia

CMFMCConvection()

GEOS-Chem RAS / Grell-Freitas convective transport operator.

No struct fields — the forcing arrays (cmfmc, optionally dtrain) arrive via TransportModel.convection_forcing, populated each substep by DrivenSimulation._refresh_forcing! from sim.window.convection.

Basis convention

The operator is basis-polymorphic. cmfmc and dtrain must be on the SAME basis as state.air_mass (dry in CATRINE usage, moist if the driver supplies moist forcing). The driver is responsible for basis correction upstream.

Fields required on ConvectionForcing

forcing.cmfmc :: AbstractArray{FT, 3} at interfaces, shape (Nx, Ny, Nz+1). Units kg / m² / s on the state's basis.
forcing.dtrain :: Union{AbstractArray{FT, 3}, Nothing} at centers, shape (Nx, Ny, Nz). When nothing, the kernel runs Tiedtke-style single-flux transport.

Face-indexed ReducedGaussian uses (ncell, Nz+1) / (ncell, Nz). CubedSphere uses NTuple{6} of per-panel (Nc, Nc, Nz+1) / (Nc, Nc, Nz) arrays.

CFL sub-cycling

The kernel sub-cycles internally based on the CMFMC profile:

julia

n_sub = max(1, ceil(max_over_grid(cmfmc × dt / bmass) / 0.5))

n_sub is cached per window in workspace.cached_n_sub[]; the driver clears the cache on window roll via invalidate_cmfmc_cache!(workspace). Bit-exact: one call with dt matches n_sub manual calls with sdt = dt / n_sub.

Well-mixed sub-cloud layer

Applies GCHP's pressure-weighted well-mixed treatment below cloud base (convection_mod.F90:742-782). Absent in legacy Julia; deliberate improvement for surface-source tracers. See the kernel Pass-0 block in cmfmc_kernels.jl.

Scope

CMFMCConvection now supports structured LatLon, face-indexed ReducedGaussian, and panel-native CubedSphere runtime state. The CS path keeps forcing panel-native too: forcing.cmfmc and forcing.dtrain are NTuple{6} payloads loaded by the CS transport driver and applied column-locally on the halo-free panel interior.

Adjoint path

The forward operator is linear in tracer mixing ratio (verified by the adjoint-identity test in test/test_cmfmc_convection.jl). NO positivity clamp is applied inside the kernel: the two-term tendency cmfmc · (q_above - q_env) + dtrain · (qc - q_env) stays linear, and tiny negativities that arise from inconsistent met data are absorbed by the global mass fixer. A future adjoint kernel reverses the two-pass order (tendency first, then updraft accumulation) with transposed coefficients; the four-term scavenging-restoring form is a wet-deposition follow-up.

AtmosTransport.Operators.Convection.CMFMCWorkspace Type

julia

CMFMCWorkspace{FT, QC, CA}

Per-sim pre-allocated workspace for CMFMCConvection.

Fields

qc_scratch :: QC — updraft concentration buffer, one entry per cell, shape matching air_mass ((Nx, Ny, Nz) for structured grids, (ncells, Nz) for face-indexed grids). Reused across all substeps and all tracers.
cell_metrics :: CA — pre-adapted cell-area metric vector used by the CFL scan and the kernel: cell_areas_by_latitude(mesh) for structured grids, per-cell cell_area(mesh, c) for face-indexed grids.
cached_n_sub :: Base.RefValue{Int} — the CFL-derived sub-step count for the current met window. Stays valid as long as cache_valid[] == true; re-computed from the current CMFMC / air-mass state on the next apply! call when invalid. Int (not Float) because it's the integer sub-step count.
cache_valid :: Base.RefValue{Bool} — sentinel; cleared via invalidate_cmfmc_cache! when the met window advances (DrivenSimulation._maybe_advance_window!, Commit 8).

Usage

Constructed once at DrivenSimulation setup time and reused for the whole run. Adapt.adapt_structure preserves the cached scalar state on the host while adapting qc_scratch to the requested backend.

AtmosTransport.Operators.Convection.NoConvection Type

julia

NoConvection()

Identity operator — apply! is a no-op. Default for configurations without active convection. Dispatch is a compile-time dead branch in TransportModel.step!, so the convection block collapses to zero floating-point work when no operator is installed (bit-exact backward- compatible with pre-plan-18 behavior).

AtmosTransport.Operators.Convection.TM5Convection Type

julia

TM5Convection(; tile_workspace_gib = 1.0)

TM5-style convective transport operator. Four-field mass-flux scheme following Tiedtke (1989) as implemented in TM5-4DVAR: two entrainment and two detrainment fields (updraft + downdraft). The backward-Euler transport matrix conv1 = I - dt·D is dense within the cloud window and identity above; the solver assembles and factorizes only the active lower-right cloud block and stores the pivot vector for adjoint replay in plan 19.

The forcing arrays (entu, detu, entd, detd) arrive via TransportModel.convection_forcing.tm5_fields, populated each substep by DrivenSimulation._refresh_forcing! from sim.window.convection.

Memory budget

tile_workspace_gib (binary GiB) is the per-topology target for the TM5 column-tile workspace. _convection_workspace_for(::TM5Convection, ...) reads this field and derives a tile size B via derive_tile_columns; the TM5Workspace then allocates a single (Nz, Nz, B) matrix slab plus matching pivot / cloud-dim / amu / amd slabs. A larger budget means fewer kernel launches per substep at the cost of larger GPU working set; the default 1.0 GiB fits all production resolutions through C720/L137 with slack on H100. The tile machinery is the load-bearing storage change: the workspace no longer scales with N_cells × Nz².

Basis convention

TM5Convection is basis-polymorphic, identical to CMFMCConvection. The four forcing fields must be on the same basis as state.air_mass (moist by upstream Fortran convention and by the ec2tm preprocessor default; dry requires a sibling preprocessor path).

Fields required on ConvectionForcing

forcing.tm5_fields :: NamedTuple{(:entu, :detu, :entd, :detd)} with all four arrays at layer centers in AtmosTransport orientation (k=1=TOA, k=Nz=surface). Units kg / m² / s. Shapes per topology:
- Structured LatLon: (Nx, Ny, Nz) per field.
- Face-indexed ReducedGaussian: (ncells, Nz) per field.
- Panel-native CubedSphere: NTuple{6, AbstractArray{FT, 3}} per field, with per-panel shape (Nc, Nc, Nz).

Orientation conversion + sign flip on entd happen in the preprocessor (src/Preprocessing/tm5_convection_conversion.jl). The operator performs zero runtime orientation gymnastics.

Solver class

Partial-pivot Gaussian elimination on the active lower-right block per column (see test/test_tm5_sparsity_above_icltop.jl for the structure survey). Identity rows above the effective cloud top and the lower-left zero quadrant are skipped by both the factorization and the tracer solve.

Pivoting is kept even though the matrix is diagonally dominant by construction (upstream Fortran comment says pivoting "not needed"). Per plan 23 principle 3, the pivot vector is stored in TM5Workspace so plan 19 (adjoint) can replay the same factorization with trans='T'.

CFL sub-cycling

None. The backward-Euler matrix solve is unconditionally stable for any dt, unlike CMFMCConvection's forward-Euler two-pass update which requires sub-cycling when the CMFMC profile is strong. The kernel launches once per tile and calls synchronize(backend) once per apply!.

AtmosTransport.Operators.Convection.TM5Workspace Type

julia

TM5Workspace{FT, M, P, C, F, A, CA}

Per-sim pre-allocated workspace for TM5Convection.

Fields

conv1 :: M — conv1 = I - dt·D matrix slab, one (Nz, Nz) block per column. Parametric on array type so Adapt.adapt_structure can swap CPU ↔ GPU without changing the TM5Workspace type constructor (plan 23 principle 5). Shapes per topology:
- Structured LatLon: (Nz, Nz, Nx, Ny) — 4D.
- Face-indexed ReducedGaussian: (Nz, Nz, ncells) — 3D.
- Panel-native CubedSphere: NTuple{6, AbstractArray{FT, 4}} with per-panel shape (Nz, Nz, Nc, Nc).
pivots :: P — permutation vector from partial-pivot LU, one Nz-length Int slice per column. Per plan 23 principle 3, preserved so plan 19 (adjoint) can replay the same factorization with transposed back-substitution. Shapes strip the leading Nz from conv1 shape: (Nz, Nx, Ny) / (Nz, ncells) / NTuple{6, (Nz, Nc, Nc)}.
cloud_dims :: C — per-column (icltop, iclbas, icllfs) triple in AtmosTransport indexing (k=1=TOA, k=Nz=surface; the preprocessor delivers forcings in this orientation so the solver has zero orientation logic). Shape (3, Nx, Ny) / (3, ncells) / NTuple{6, (3, Nc, Nc)}.
f_scratch :: F — per-column intermediate matrix for the matrix build (TM5 f(0:lmx, 1:lmx) after the storage plan's Commit 3 merged the updraft into f). Plan 23 Commit 4 pre-allocates this so _tm5_build_conv1! runs without heap allocation inside KA kernels (mandatory on GPU; same contract on CPU for parity). f_scratch aliases conv1 in the production workspace because conv1 is only needed after f has been converted into I - dt*D; this saves one dense (Nz, Nz) slab per column. The standalone fu_scratch field that previously carried the updraft contribution was dropped in Commit 3 — the updraft and downdraft passes write disjoint index ranges, so the builder writes directly into f.
amu_scratch :: A, amd_scratch :: A — length-(Nz+1) per-column boundary-aware mass-flux vectors (TM5 amu(0:lmx) / amd(0:lmx)). Same allocation policy as the scratch matrices.
cell_metrics :: CA — pre-adapted cell-area metrics used to convert runtime air_mass from kg per grid cell to the kg/m² column mass expected by the TM5 matrix. nothing is allowed only for standalone unit-area column tests; topology-level apply! requires this field to be populated.

Usage

Constructed once at DrivenSimulation setup time via _convection_workspace_for(::TM5Convection, state, grid). Adapt.adapt_structure adapts each array to the requested backend without reallocating on the host.

AtmosTransport.Operators.Convection.TM5Workspace Method

julia

TM5Workspace(air_mass; tile_columns = …,
                        tile_workspace_gib = nothing,
                        cell_metrics = nothing) -> TM5Workspace

Construct a fresh workspace from an air-mass payload. air_mass may be a single array (structured (Nx, Ny, Nz) or face-indexed (ncells, Nz)) or a cubed-sphere panel tuple — the workspace itself is topology-agnostic. The per-column element type of conv1 matches eltype(air_mass).

The per-launch column count B is set by exactly one of:

tile_columns::Integer (explicit). Default _tm5_total_cells_per_launch(air_mass) — one tile covers the whole launch and the workspace is bit-equal to the pre-Commit-4 per-cell allocator. Production paths use this branch when they already know B.
tile_workspace_gib::Real (budget). Picks B via derive_tile_columns from TM5Convection's tile_workspace_gib field. Bypassed if tile_columns is also passed (explicit wins).
cell_metrics carries topology cell areas. Production _convection_workspace_for(::TM5Convection, ...) supplies this; leaving it nothing is only for unit-area solver tests.

Specifying both is an error.

AtmosTransport.Operators.Convection.apply_convection! Method

julia

apply_convection!(q_raw, air_mass, forcing::ConvectionForcing,
                   ::NoConvection, dt, workspace, grid) -> nothing

Array-level no-op, parallels the plan 16b / 17 pattern (apply_vertical_diffusion!, apply_surface_flux!). Accepts any shape of q_raw / air_mass — NoConvection doesn't inspect them. Returns nothing.

The plan-18 structured apply! flow goes through the state-level method above. apply_convection! is reserved for the future case where the convection block is called from inside a palindrome or another composed setting — same signature contract as plan 16b and 17.

AtmosTransport.Operators.Convection.apply_convection! Method

julia

apply_convection!(q_raw, air_mass, forcing::ConvectionForcing,
                   op::TM5Convection, dt, workspace::TM5Workspace,
                   grid::AtmosGrid) -> nothing

Array-level entry point — parallels the CMFMC contract at operators.jl:70-89. Dispatches on grid mesh type and launches the matching KA kernel from tm5_kernels.jl. Single synchronize(backend) at the end (TM5 matrix solve is unconditionally stable; no sub-cycling).

AtmosTransport.Operators.Convection.derive_tile_columns Method

julia

derive_tile_columns(::Type{FT}, Nz, budget_gib, total_cells) -> Int

Pick a tile size B from a per-cell memory cost and a target budget. The cost model accounts for every per-cell field that TM5Workspace tiles:

conv1 : Nz² × sizeof(FT)
amu_scratch + amd_scratch : 2(Nz+1) × sizeof(FT)
pivots : Nz × sizeof(Int)
cloud_dims : 3 × sizeof(Int)

f_scratch is a structural alias for conv1 (saves one matrix slab per cell). B is clamped between 256 (avoid pathological launches) and total_cells (one tile covers the whole topology when the budget allows it).

AtmosTransport.Operators.Convection.invalidate_cmfmc_cache! Method

julia

invalidate_cmfmc_cache!(ws::CMFMCWorkspace) -> nothing

Mark the cached CFL n_sub as stale. Called on met-window advance by DrivenSimulation._maybe_advance_window! (plan 18 Commit 8). The next apply! recomputes n_sub from the fresh CMFMC array.

AtmosTransport.Operators.Convection.invalidate_tm5_cache! Method

julia

invalidate_tm5_cache!(ws::TM5Workspace) -> nothing
invalidate_tm5_cache!(::Any) -> nothing

Mark the P6 LU-factor cache stale. Called on met-window advance by DrivenSimulation._maybe_advance_window! so the next apply! repopulates the cache from the new forcing fields. The polymorphic no-op default keeps the call site type-agnostic (the same way invalidate_cmfmc_cache! is structured).

When the cache is disabled (ws.cache_valid === nothing) this is also a no-op — the workspace simply has nothing to invalidate.

Diffusion

AtmosTransport.Operators.Diffusion.AbstractSurfaceFluxCoupling Type

julia

AbstractSurfaceFluxCoupling

Typed policy for where configured surface fluxes enter relative to a vertical diffusion/mixing operator.

SplitSurfaceFluxCoupling: existing Strang-center composition, V(dt/2) -> S(dt) -> V(dt/2).
DiffusiveSurfaceFluxBoundary: GCHP/VDIFF-style lower-boundary placement, S(dt) -> V(dt), so fresh surface flux is included in the implicit vertical mixing solve.

AtmosTransport.Operators.Diffusion.ImplicitVerticalDiffusion Type

julia

ImplicitVerticalDiffusion(; kz_field)

Backward-Euler vertical diffusion driven by a cell-centered Kz field. Two spatial layouts are supported:

structured: AbstractTimeVaryingField{FT, 3} over (Nx, Ny, Nz)
face-indexed: AbstractTimeVaryingField{FT, 2} over (ncells, Nz)

Concrete examples:

ConstantField{FT, 3} / ConstantField{FT, 2}
ProfileKzField{FT} with default rank 3 or ProfileKzField(profile; spatial_rank = 2)
PreComputedKzField{FT, A} wrapping 3D or 2D storage
DerivedKzField for meteorology-driven Beljaars-Viterbo on structured grids
CubedSphereField wrapping one structured rank-3 field per panel

apply! contract

julia

apply!(state, meteo, grid, op::ImplicitVerticalDiffusion, dt; workspace)

Refreshes the Kz cache: update_field!(op.kz_field, t) with t drawn from the meteorology where available; plan 16b currently passes zero(FT) as a placeholder (chemistry-style, mirrors plan 15's deferred current_time(meteo) accessor).
Reads workspace.dz_scratch as the current layer thicknesses [m]. The caller is responsible for filling this array before calling apply! — typically from a hydrostatic integration of the current delp and surface temperature.
Uses workspace.w_scratch as Thomas-forward-elimination storage.
Launches a layout-specific diffusion kernel:
- structured: _vertical_diffusion_kernel! over (Nx, Ny, Nt)
- face-indexed: _vertical_diffusion_face_kernel! over (ncells, Nt)

The operator is linear (Kz does not depend on tracer values), so a single apply!(dt) at the palindrome center is equivalent to two half-steps — see plan 16b §4.3 Decision 8. Commit 4 performs the palindrome integration.

Fields

kz_field::KzF — any AbstractTimeVaryingField{FT, 2} or AbstractTimeVaryingField{FT, 3} providing cell-centered Kz values [m²/s geometric].

AtmosTransport.Operators.Diffusion.NoDiffusion Type

julia

NoDiffusion()

Identity operator — apply! is a no-op. Default for configurations without active vertical mixing, and the value strang_split_mt! sees when the palindrome's V position is unoccupied (Commit 4).

AtmosTransport.Operators.Diffusion._vertical_diffusion_kernel! Method

julia

_vertical_diffusion_kernel!(q, kz_field, dz, w_scratch, dt, Nz)

KernelAbstractions kernel: implicit (Backward-Euler) vertical diffusion for one column per (i, j, t) thread.

q::AbstractArray{FT, 4} — tracer values (Nx, Ny, Nz, Nt), read for old values and written with new values in place.
kz_field::AbstractTimeVaryingField{FT, 3} — Kz at cell centers.
dz::AbstractArray{FT, 3} — layer thicknesses in meters, (Nx, Ny, Nz). Caller supplies; not mutated.
w_scratch::AbstractArray{FT, 3} — caller-supplied workspace, (Nx, Ny, Nz). Holds the Thomas forward-elimination factors between the forward and back-substitution loops.
dt::FT — time step.
Nz::Int — number of vertical levels (passed explicitly so the kernel avoids a size call).

The (a, b, c, d) tridiagonal entries are named local values at each level k rather than pre-built arrays. The coefficient formulas exactly match build_diffusion_coefficients — the reference used by the test suite. Any change to the formulas here requires updating the reference (or vice versa).

Adjoint note

The future adjoint kernel is structurally identical: read Kz and dz the same way, build the same (a_k, b_k, c_k) per level, then apply the transposition rule at the tridiagonal interface (a_T[k] = c[k-1], b_T[k] = b[k], c_T[k] = a[k+1]) before passing to a Thomas solve. Archival adjoint template: docs/resources/developer_notes/legacy_adjoint_templates/boundary_layer_diffusion_adjoint.jl:74-84.

AtmosTransport.Operators.Diffusion.apply_vertical_diffusion! Function

julia

apply_vertical_diffusion!(q_raw, op, workspace, dt, meteo = nothing) -> nothing

Low-level entry point. Applies one Backward-Euler diffusion step to a raw tracer buffer in any of the supported layouts:

structured packed tracers: q_raw :: (Nx, Ny, Nz, Nt)
face-indexed packed tracers: q_raw :: (ncells, Nz, Nt)
face-indexed single-tracer slice: q_raw :: (ncells, Nz)

This is the function strang_split_mt! calls at the palindrome center (plan 16b Commit 4). The face-indexed reduced-Gaussian path also uses it at its H → V → D → V → H center slot.

meteo is threaded through to update_field!(op.kz_field, t) as t = FT(current_time(meteo)) (or zero(FT) if meteo === nothing). Plan 17 Commit 4: meteo defaults to nothing so pre-17 palindrome call sites (apply_vertical_diffusion!(rm, op, ws, dt)) continue to work unchanged; Commit 5 threads meteo through the palindrome.

NoDiffusion is a no-op: the method is = nothing so Julia's dispatch reduces the call site to a dead branch when diffusion_op isa NoDiffusion. This is what makes the Commit 4 palindrome integration bit-exact backward-compatible.

AtmosTransport.Operators.Diffusion.apply_vertical_diffusion_vmr! Method

julia

apply_vertical_diffusion_vmr!(rm, air_mass, op, workspace, dt, meteo; halo_width)

Cubed-sphere helper for state variables stored as tracer mass. The implicit vertical solver acts on mixing ratio; this wrapper converts tracer mass to VMR using the current dry air mass, applies the existing column solve, then restores tracer mass before advection resumes.

AtmosTransport.Operators.Diffusion.build_diffusion_coefficients Method

julia

build_diffusion_coefficients(Kz_col, dz_col, dt) -> (a, b, c)

Reference Backward-Euler tridiagonal coefficient builder for vertical diffusion on one column. Returns three Nz-length Vector{FT}s representing T = I - dt · D, where D is the discrete second-derivative operator with variable Kz at cell centers and Neumann (zero-flux) boundary conditions.

Interface Kz is the arithmetic mean of adjacent cell-center values; interface dz is the arithmetic mean of adjacent cell thicknesses. Zero-flux BCs set a[1] = 0 and c[Nz] = 0.

Formulas at interior k:

julia

dz_above = 0.5 × (dz[k-1] + dz[k])
dz_below = 0.5 × (dz[k] + dz[k+1])
Kz_above = 0.5 × (Kz[k-1] + Kz[k])
Kz_below = 0.5 × (Kz[k] + Kz[k+1])
D_above  = Kz_above / (dz[k] × dz_above)
D_below  = Kz_below / (dz[k] × dz_below)
a[k]     = -dt × D_above
b[k]     =  1 + dt × (D_above + D_below)
c[k]     = -dt × D_below

At the boundaries: D_above = 0 at k = 1, D_below = 0 at k = Nz.

Role: reference vs. production

This function is the reference used in tests. The production kernel _vertical_diffusion_kernel! inlines the same formulas at each level (to avoid allocation and to read Kz through field_value). Tests verify the kernel's output matches the output of this reference on matched inputs.

If the formulas here are ever changed, the kernel must be kept in lock-step.

Adjoint note

As in solve_tridiagonal!, the transpose of the resulting tridiagonal is a_T[k] = c[k-1], b_T[k] = b[k], c_T[k] = a[k+1]. The formulas above are symmetric in Kz / dz (the interface values), so the adjoint operator built by a future plan differs only by the shift, not by the coefficient content.

AtmosTransport.Operators.Diffusion.fill_dz_hydrostatic_constT! Method

julia

fill_dz_hydrostatic_constT!(dz, ps, ak_ifc, bk_ifc;
                             T_ref = 260, R = 287.04, gravity = 9.81)

Populate a 3D (Nx, Ny, Nz) dz array (host or device) from surface pressure ps and hybrid sigma-pressure interface coefficients ak_ifc, bk_ifc (length Nz + 1). Backend follows dz.

AtmosTransport.Operators.Diffusion.fill_dz_hydrostatic_constT! Method

julia

fill_dz_hydrostatic_constT!(dz::AbstractArray{<:Any, 2},
                             ps::AbstractArray{<:Any, 1},
                             ak_ifc, bk_ifc; ...)

Face-indexed variant for ReducedGaussian topology: dz shape (ncells, Nz), ps shape (ncells,). Same constant-T_ref formula as the structured/CS overloads, just unrolled over the face-indexed cell axis.

AtmosTransport.Operators.Diffusion.fill_dz_hydrostatic_constT! Method

julia

fill_dz_hydrostatic_constT!(dz_panels::NTuple{6}, ps_panels::NTuple{6},
                             ak_ifc, bk_ifc; ...)

Cubed-sphere variant: per-panel 3D (Nc, Nc, Nz) dz arrays are filled from per-panel ps arrays of shape (Nc, Nc) (interior only — the runtime stores surface_pressure without the advection halo).

AtmosTransport.Operators.Diffusion.solve_tridiagonal! Method

julia

solve_tridiagonal!(x, a, b, c, d, w)

Solve the tridiagonal linear system T x = d in place, where T has sub-diagonal a, main diagonal b, super-diagonal c. By convention a[1] and c[Nz] are ignored (no-neighbor positions).

Argument roles:

x::AbstractVector{FT} — output, overwritten with the solution.
a, b, c, d::AbstractVector{FT} — read only, not mutated.
w::AbstractVector{FT} — caller-supplied workspace, length ≥ Nz. Holds the Thomas forward-elimination factors w[k] = c[k] / denom used during back-substitution.

Implements the standard Thomas algorithm with one extra array (w) so that b and d can be read-only. For per-column Nz this is Θ(Nz) arithmetic, no allocation.

Adjoint note (forward-only in this commit)

The matrix transpose T^T is obtained by swapping shifted sub- and super-diagonals:

julia

a_T[k] = c[k - 1]     # for k ≥ 2
b_T[k] = b[k]
c_T[k] = a[k + 1]     # for k ≤ Nz - 1

A future adjoint kernel calls this same solve_tridiagonal! after building (a_T, b_T, c_T). No structural change required. Archival adjoint template: docs/resources/developer_notes/legacy_adjoint_templates/boundary_layer_diffusion_adjoint.jl:74-84.

SurfaceFlux

AtmosTransport.Operators.SurfaceFlux.AbstractSurfaceFluxOperator Type

julia

AbstractSurfaceFluxOperator

Top of the surface emission operator hierarchy. Concrete subtypes in plan 17 Commit 3:

NoSurfaceFlux — identity (default).
SurfaceFluxOperator — wraps a PerTracerFluxMap and applies each tracer's surface flux to k = Nz of the tracer mass array.

Every concrete subtype implements two entry points:

State-level apply!(state, meteo, grid, op, dt; workspace) → mutates state.tracers_raw.
Array-level apply_surface_flux!(q_raw, op, ws, dt, meteo, grid; tracer_names) → mutates a raw tracer buffer directly. Supported layouts are:
structured packed: (Nx, Ny, Nz, Nt)
face-indexed packed: (ncells, Nz, Nt)
face-indexed single-tracer slice: (ncells, Nz)
cubed-sphere single-tracer panels: NTuple{6} of (Nc + 2Hp, Nc + 2Hp, Nz) Used by the structured multi-tracer palindrome (Commit 5) and the reduced-Gaussian face-indexed transport block (plan 22A).

AtmosTransport.Operators.SurfaceFlux.NoSurfaceFlux Type

julia

NoSurfaceFlux()

Identity operator — apply! is a no-op. Default for configurations without surface emissions, and the value strang_split_mt! sees when the palindrome's S position is unoccupied (Commit 5).

NoSurfaceFlux's apply! is literally = state (and the array-level apply_surface_flux! is = nothing). Julia's dispatch turns the call site into a dead branch — zero floating-point work, bit-exact backward compatibility.

AtmosTransport.Operators.SurfaceFlux.PerTracerFluxMap Type

julia

PerTracerFluxMap{S <: Tuple}

An ordered tuple of SurfaceFluxSources, keyed by the tracer_name field on each entry. Supplies the surface flux operator (Commit 3) with per-tracer surface-rate data and ensures efficient tuple-splatting at kernel-launch time.

The map is NTuple-backed rather than Dict-backed. Rationale:

Existing DrivenSimulation.surface_sources::Tuple storage is already a tuple; the map type is a thin re-export with a lookup helper.
Tuples are bits-stable on GPU (captured by kernels without boxing / hashing); Dicts would require special Adapt handling and per-launch lookup cost.
Small-N (typically 1-5 emitting tracers) makes linear scan cheaper than hashing anyway.

Tracers absent from the map have zero surface flux. The consumer operator iterates the map and applies each source; it does NOT walk the state's full tracer_names and look up each by name.

Construction

julia

co2_rate  = fill(2.0e-7, Nx, Ny)    # kg/s per cell for :CO2
sf6_rate  = fill(1.5e-9, Nx, Ny)    # kg/s per cell for :SF6
rn222_rate = fill(3.0e-11, Nx, Ny)  # kg/s per cell for :Rn222

map = PerTracerFluxMap(
    SurfaceFluxSource(:CO2,   co2_rate),
    SurfaceFluxSource(:SF6,   sf6_rate),
    SurfaceFluxSource(:Rn222, rn222_rate),
)

length(map)                 # 3
flux_for(map, :CO2) === ... # the :CO2 SurfaceFluxSource
flux_for(map, :CH4)         # nothing (not emitting)

Adapt / GPU

Adapt.adapt(CuArray, map) converts each cell_mass_rate to a CuArray transparently via SurfaceFluxSource's adapt_structure; tuples adapt element-wise.

Fields

sources :: S — NTuple{N, SurfaceFluxSource{<:AbstractArray}} for some N. Each entry carries a tracer name and the rate array.

AtmosTransport.Operators.SurfaceFlux.PerTracerFluxMap Method

julia

PerTracerFluxMap(sources::AbstractVector{<:SurfaceFluxSource})
PerTracerFluxMap(sources::Tuple)

Generic collection constructor. The input is frozen into an NTuple.

AtmosTransport.Operators.SurfaceFlux.PerTracerFluxMap Method

julia

PerTracerFluxMap(sources::SurfaceFluxSource...)

Variadic constructor: PerTracerFluxMap(src1, src2, src3) wraps the three sources into an NTuple-backed map.

AtmosTransport.Operators.SurfaceFlux.SurfaceFluxOperator Type

julia

SurfaceFluxOperator(flux_map::PerTracerFluxMap)
SurfaceFluxOperator(sources::SurfaceFluxSource...)

Applies a PerTracerFluxMap of surface sources to the k = Nz slab of the tracer mass array.

For every tracer named in the flux map, the operator launches the layout-appropriate surface-flux kernel and adds rate × dt to the surface layer k = Nz of the matching tracer:

structured packed state (Nx, Ny, Nz, Nt) → _surface_flux_kernel! over (Nx, Ny)
face-indexed packed state (ncells, Nz, Nt) → _surface_flux_face_kernel! over ncells

Tracer indices are resolved on the host from state.tracer_names. Tracers absent from the map are untouched.

Per plan 17 Decision 1, the rate is in kg/s per cell (already area- integrated); no cell-area multiplier appears in the kernel. This preserves the pre-17 _apply_surface_source! semantics and matches the way legacy SurfaceFluxSource callers supply their arrays.

apply! contract

julia

apply!(state, meteo, grid, op::SurfaceFluxOperator, dt; workspace)

Walks op.flux_map in storage order.
For each source, skips tracers not present in state.
Dispatches to the array-level entry point with tracer_names pulled from the state.
Synchronises the backend once at the end.

Fields

flux_map :: M — a PerTracerFluxMap of SurfaceFluxSources. Emitting tracers are exactly those named in the map.

AtmosTransport.Operators.SurfaceFlux.SurfaceFluxOperator Method

julia

SurfaceFluxOperator(sources::SurfaceFluxSource...)

Convenience constructor: wraps a variadic list of sources in a PerTracerFluxMap first. Empty variadic list is allowed (produces an empty-map operator; equivalent to NoSurfaceFlux except type- distinguishable).

AtmosTransport.Operators.SurfaceFlux.SurfaceFluxSource Type

julia

SurfaceFluxSource{RateT}

A single-tracer surface source: a tracer_name plus a cell_mass_rate array supplying mass added per cell per second to the surface layer.

tracer_name :: Symbol — matches a name in CellState.tracer_names.
cell_mass_rate :: RateT — one of:
- a 2D (Nx, Ny) array for structured grids
- a 1D (Nc,) array for face-indexed grids
- an NTuple{6} of 2D (Nc, Nc) arrays for cubed-sphere panels
The units are kg/s per cell — already area-integrated. The surface flux kernel applies rm_surface += rate × dt without multiplying by cell area.

Why kg/s per cell (not kg/m²/s)

Plan 17 Decision 1 (user-approved 2026-04-19): the prognostic tracer is mass (kg per cell), so a per-cell rate × dt is the natural unit and matches the legacy DrivenSimulation._apply_surface_source! signature that shipped before plan 17. A per-area (kg/m²/s) variant that multiplies by cell_area is deferred to a follow-up plan.

History

Originally introduced in src/Models/DrivenSimulation.jl before plan 17. Migrated to src/Operators/SurfaceFlux/ in plan 17 Commit 2 so that SurfaceFluxOperator (Commit 3) can consume it; the name is still re-exported from AtmosTransport for backward compat with external callers that imported it by fully-qualified name.

Fields

tracer_name :: Symbol
cell_mass_rate :: RateT — backend-agnostic; Adapt.adapt converts the array between host and device transparently via Adapt.adapt_structure.

AtmosTransport.Operators.SurfaceFlux.apply_surface_flux! Function

apply_surface_flux!(q_raw, op::AbstractSurfaceFluxOperator, workspace, dt, meteo, grid; tracer_names)

Array-level surface-flux application. Writes directly to the supplied tracer buffer q_raw, adding each source's rate × dt contribution to the surface slab k = Nz.

Supported layouts:

structured packed: q_raw :: (Nx, Ny, Nz, Nt)
face-indexed packed: q_raw :: (ncells, Nz, Nt)
face-indexed single-tracer slice: q_raw :: (ncells, Nz)

tracer_names::NTuple{Nt, Symbol} is required as a keyword so the function can resolve each source's name to a slab index without reaching back into the caller's CellState. This lets the palindrome integration (Commit 5) point the operator at either the caller's state.tracers_raw or the workspace's ping-pong buffer — whichever currently holds the post-Z-sweep tracer state.

workspace, meteo, and grid are accepted but currently unused; they are in the signature to match the operator-interface convention (§"Workflow: Adding a new physics operator" in CLAUDE.md) and to leave room for future extensions (e.g. meteorology-dependent emissions).

Returns nothing on success. For NoSurfaceFlux, returns nothing immediately (zero floating-point work).

AtmosTransport.Operators.SurfaceFlux.emitting_tracer_indices Method

julia

emitting_tracer_indices(op::SurfaceFluxOperator, state::CellState) -> NTuple

Ordered tuple of tracer indices in state.tracer_names for the tracers present in op.flux_map. Tracers in the map but missing from the state are skipped (returned as nothing slots). Useful for testing / introspection and unchanged by apply!.

AtmosTransport.Operators.SurfaceFlux.flux_for Method

julia

flux_for(map, tracer_name::Symbol) -> SurfaceFluxSource | nothing

Return the SurfaceFluxSource for the named tracer, or nothing if the tracer has no surface source in this map. O(N) linear scan, which is fine for typical N ≤ 10.

Chemistry

AtmosTransport.Operators.Chemistry.CompositeChemistry Type

julia

CompositeChemistry(schemes...)
CompositeChemistry(schemes::Tuple)

Apply multiple chemistry operators sequentially. Used when different species need independent transformations or when different operator types (decay + photolysis + ...) must run in a prescribed order.

julia

chem = CompositeChemistry(
    ExponentialDecay(; Rn222 = 330_350.4),
    ExponentialDecay(; Kr85  = 3.394e8),
)

AtmosTransport.Operators.Chemistry.ExponentialDecay Type

Keyword constructor: ExponentialDecay(; Rn222 = half_life_seconds, ...).

AtmosTransport.Operators.Chemistry.ExponentialDecay Type

julia

ExponentialDecay{FT, N, R}(decay_rates, tracer_names)

Multi-tracer first-order decay: c *= exp(-rate * dt) applied in-place to every selected tracer at every cell. Exact for constant rate and any dt; unconditionally stable; trivially parallel.

Fields

decay_rates :: R — an NTuple{N, <: AbstractTimeVaryingField{FT, 0}} of rate-valued fields, one per selected tracer [1/s]. apply! calls update_field! on each rate before launching the kernel.
tracer_names :: NTuple{N, Symbol} — which tracers this operator applies to

Construction

julia

ExponentialDecay(; Rn222 = 330_350.4)                   # from half-lives [s]
ExponentialDecay(Float32; Rn222 = 330_350.4, Kr85 = 3.394e8)

The keyword constructor converts half-life T to decay rate λ = log(2) / T (first-order exponential decay) and wraps each rate in a ConstantField{FT, 0}. Future plans may pass time-varying rates (e.g. temperature-dependent reaction rates) through the same field interface.

Common isotopes:

²²²Rn: half-life = 330_350.4 s (3.8235 days) → λ ≈ 2.098e-6 s⁻¹
⁸⁵Kr: half-life = 3.394e8 s (10.76 years) → λ ≈ 2.042e-9 s⁻¹

AtmosTransport.Operators.Chemistry.NoChemistry Type

julia

NoChemistry()

Identity operator — apply! is a no-op. Default for runs without active chemistry.

AtmosTransport.Operators.Chemistry.chemistry_block! Function

julia

chemistry_block!(state::CellState, meteo, grid, operators, dt;
                 workspace=nothing)

Apply each operator in the operators tuple to state, in order. Returns state for chaining.

The operators collection may be:

a Tuple of AbstractChemistryOperators (zero or more),
a single AbstractChemistryOperator (wrapped into a 1-tuple internally),
an empty Tuple() (identity).

This is the step-level chemistry entry point; apply!(::CompositeChemistry) is the operator-level composer. Both exist and both are supported.