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 Module

julia

Operators

Physics operators for the basis-explicit transport architecture.

Provides:

Abstract operator types (AbstractOperator, AbstractDiffusion, etc.)
Advection hierarchy: AbstractAdvectionScheme → AbstractConstantScheme, AbstractLinearScheme, AbstractQuadraticScheme with concrete UpwindScheme, SlopesScheme, and structured-grid PPMScheme
Strang splitting orchestrator: strang_split!, apply!

source

AtmosTransport.Operators.AbstractConvection Type

julia

AbstractConvection <: AbstractOperator

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

source

AtmosTransport.Operators.AbstractDiffusion Type

julia

AbstractDiffusion <: AbstractOperator

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

source

AtmosTransport.Operators.AbstractOperator Type

julia

AbstractOperator

Root type for all physics operators in the transport model.

source

Advection

AtmosTransport.Operators.Advection Module

julia

Advection

Advection operators for the basis-explicit transport architecture.

Provides:

Scheme hierarchy:

UpwindScheme <: AbstractConstantScheme — first-order upwind via generic kernels
SlopesScheme <: AbstractLinearScheme — van Leer slopes (limiter-dispatched)
PPMScheme <: AbstractQuadraticScheme — structured-grid PPM (not yet an official real-data reference path)
LinRoodPPMScheme <: AbstractAdvectionScheme — cubed-sphere FV3/Lin-Rood PPM path
AbstractLimiter subtypes: NoLimiter, MonotoneLimiter, PositivityLimiter

Multi-tracer optimization:

TracerView — zero-cost 3D slice adapter for 4D tracer arrays
Multi-tracer kernel shells fuse the N-tracer loop into GPU kernels, reducing launches from 6N to 6 per Strang split

Infrastructure:

AdvectionWorkspace + strang_split! — Strang splitting orchestrator
CFL utilities for subcycling decisions

source

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)

source

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

source

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.

source

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

source

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

source

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.

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

source

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.

source

AtmosTransport.Operators.Advection.CSAdvectionWorkspace Type

julia

CSAdvectionWorkspace{FT, A, P3, A4, P4}

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.
m_pp_buf, rm_4d_pp_buf are full-panel spare buffers for the packed ping-pong path, avoiding the per-sweep copy-back kernels.
max_subcycles tracks this workspace's high-water mark for CFL diagnostics; keeping it with the workspace prevents unrelated simulations sharing state.

source

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

source

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.

source

AtmosTransport.Operators.Advection.NoAdvection Type

julia

NoAdvection <: AbstractAdvectionScheme

Identity scheme — apply! is a no-op. Default for runs that need to isolate other operators (convection, chemistry, diagnostics) from the advective transport block, e.g. a "convection-alone" timing experiment or a synthetic-state regression.

When NoAdvection is selected:

The advection sweeps are skipped on all three topologies (LatLonMesh, ReducedGaussianMesh, CubedSphereMesh).
Diffusion is allowed: with no advection to wrap, the V(dt) diffusion step is applied directly through the mass-flux VMR kernels, which are mass-conserving on their own and need no Strang palindrome. This is the natural "diffusion-only" experimental setup.
Surface emissions are rejected with an actionable ArgumentError. They are integrated as Strang half-steps wrapping the advection block, so running them without advection silently drops 2nd-order accuracy. To run a convection-only configuration omit [tracers.*.surface_flux] blocks.
The CS workspace is auto-allocated by the TransportModel constructor when diffusion is active (and is nothing when diffusion is also off).

Example

julia

scheme = NoAdvection()

source

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.

source

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)

source

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

source

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

source

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]

source

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

source

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.

source

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.

source

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.

source

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.

source

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.

source

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.

source

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.

source

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.

source

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

source

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.

source

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

source

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

source

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.

source

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.

source

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

source

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

source

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

source

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

source

Convection

AtmosTransport.Operators.Convection Module

julia

Convection

Convective transport operators.

This submodule ships the convection-operator hierarchy:

AbstractConvection — supertype for all convective transport operators.
NoConvection — identity no-op; default for configurations without active convection. Dispatch is a compile-time dead branch in TransportModel.step! so the presence of the convection block has zero cost for users who don't opt in.

Operators:

ConvectionForcing in ..MetDrivers carries the convective forcing arrays (with copy_convection_forcing!, allocate_convection_forcing_like, and window-struct integration).
CMFMCConvection (GCHP path, CMFMC+DTRAIN kernel with mandatory CFL sub-cycling and well-mixed sub-cloud).
TM5Convection (TM5 column solver) with TM5Workspace, running on all three topologies (LL / RG / CS).

The step!-level runtime block wires convection across all three topologies (LL / RG / CS).

apply! contract

julia

apply!(state::CellState{B},
       forcing::ConvectionForcing,
       grid::AtmosGrid,
       op::AbstractConvection,
       dt::Real;
       workspace) where {B <: AbstractMassBasis}

The operator takes ConvectionForcing directly (not a transport window or driver). _refresh_forcing! populates model.convection_forcing each substep by copying from sim.window.convection. TransportModel.step! executes the convection block between transport and chemistry. No meteo kwarg — the forcing arrays are the time information; the operator does not call current_time.

Face-indexed scope

CMFMCConvection now supports structured LatLon, face-indexed ReducedGaussian, and panel-native CubedSphere state layouts. The CS path keeps forcing panel-native too: the driver loads cmfmc / dtrain as per-panel tuples and the operator applies the same column-local logic on each panel interior.

TM5Convection ships the same three-topology scope.

The no-op NoConvection path accepts any state shape — it's a pure dead branch.

source

AtmosTransport.Operators.Convection.CMFMCConvection Type

julia

CMFMCConvection(; clamp=false)

Explicit, conservative convection driven by per-layer convective mass flux and detrainment from ConvectionForcing. Supported on lat-lon, reduced-Gaussian, and cubed-sphere states. Setting clamp=true applies the positivity correction and conservative column rescaling used for strong-CFL forcing.

source

AtmosTransport.Operators.Convection.CMFMCMatrixConvection Type

julia

CMFMCMatrixConvection(; tile_workspace_gib = 1.0,
                        use_collab_lu = false,
                        lmax_conv = 0,
                        n_merge = 1)

GEOS-derived rates through the TM5 LU solver. The kwarg-tuned knobs all pass through to the inner TM5Convection that performs the column solve — see its docstring for tile_workspace_gib, use_collab_lu, lmax_conv, and n_merge semantics.

Forcing requirements

forcing.cmfmc and forcing.dtrain must be populated on the GEOS binary's basis (whatever state.air_mass is — dry for our production runtime). forcing.tm5_fields is ignored. The operator runs once-per-met-window derivation to convert (cmfmc, dtrain) → (entu, detu, 0, 0), caching the result in CMFMCMatrixWorkspace for reuse across all substeps.

Conservation contract

Σ(m·q) is preserved to floating-point roundoff for any inert tracer: the TM5 LU matrix is column-stochastic once Σ entu = Σ detu, which the derivation kernels guarantee by absorbing any boundary-residual at the surface layer.

Adjoint

The derivation is constant in the state vector — it depends only on forcings — so the state-space adjoint is exactly the inner TM5 LU adjoint applied to the cached derived rates. See Adjoints/ConvectionAdjoint.jl.

source

AtmosTransport.Operators.Convection.CMFMCMatrixWorkspace Type

julia

CMFMCMatrixWorkspace{FT, TS, EU, ED, V}

Per-sim pre-allocated workspace for CMFMCMatrixConvection.

Fields

tm5_workspace :: TS — the underlying TM5Workspace carrying the column-tile matrix slab, pivots, and TM5 cell-area metrics.
derived_entu :: EU — derived updraft entrainment rate, same shape and basis as forcing.dtrain (LL: (Nx,Ny,Nz), RG: (ncells,Nz), CS: NTuple{6, (Nc,Nc,Nz)}). Refreshed on every met-window advance via invalidate_cmfmc_matrix_cache!; reused bit-exact across all substeps within a window.
derived_detu :: EU — derived updraft detrainment rate, same shape as derived_entu. Includes the explicit GEOS dtrain plus any negative-entrainment residual folded in to keep entu ≥ 0.
zero_entd :: ED — immutable zeros, same shape as derived_entu. Passed as the TM5 downdraft entrainment field (CMFMC has no downdraft in our binary contract).
zero_detd :: ED — immutable zeros, same shape (TM5 downdraft detrainment).
derived_valid :: V — Base.RefValue{Bool} flagging whether the cached derivation matches the current met window. Cleared via invalidate_cmfmc_matrix_cache! on window advance.

Lifecycle

The derived rates are computed once per met-window from forcing.cmfmc + forcing.dtrain by _launch_cmfmc_matrix_derivation! (see cmfmc_matrix_kernels.jl). Within a window, apply! reuses the cached arrays directly — the LU solve sees the same (entu, detu, 0, 0) quadruplet every substep.

source

AtmosTransport.Operators.Convection.CMFMCMatrixWorkspace Method

julia

CMFMCMatrixWorkspace(air_mass; tile_workspace_gib = nothing,
                      tile_columns = nothing,
                      cell_metrics = nothing,
                      halo_width = 0) -> CMFMCMatrixWorkspace

Construct a fresh workspace from an air-mass payload. air_mass may be a single array (LL / RG) or a panel tuple (CS). The inner TM5Workspace is built with the same tile-budget rules as a bare TM5Workspace; the derived-rate slabs are sized to mirror the layout of forcing.dtrain (LL/RG: same as air_mass; CS: halo-stripped from air_mass).

cell_metrics is the topology cell-area metric (mandatory in production — the TM5 matrix needs it to convert kg-per-cell air mass to the kg/m² basis the convective fluxes share).

halo_width is the per-panel halo half-width Hp for cubed-sphere payloads (the production CS factory passes grid.horizontal.Hp). LL and RG leave it at 0.

source

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

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.

source

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 the no-op behavior).

source

AtmosTransport.Operators.Convection.TM5Convection Type

julia

TM5Convection(; tile_workspace_gib=1.0, use_collab_lu=false,
                lmax_conv=0, n_merge=1)

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.

use_collab_lu selects the workgroup-collaborative solver. lmax_conv=0 uses all vertical levels; a positive value limits the active lower atmosphere. n_merge optionally aggregates adjacent active layers before the solve.

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"). The pivot vector is stored in TM5Workspace so the 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!.

source

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. 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. Preserved so the 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), with the updraft merged into f). Pre-allocated 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. There is no standalone fu_scratch field — 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.

source

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

source

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 the diffusion and surface-flux apply! pattern. Accepts any shape of q_raw / air_mass — NoConvection doesn't inspect them. Returns nothing.

The 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 the diffusion and surface-flux array entry points.

source

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

source

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

source

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!. The next apply! recomputes n_sub from the fresh CMFMC array.

source

AtmosTransport.Operators.Convection.invalidate_cmfmc_matrix_cache! Method

julia

invalidate_cmfmc_matrix_cache!(ws::CMFMCMatrixWorkspace) -> nothing
invalidate_cmfmc_matrix_cache!(::Any) -> nothing

Mark the derived (entu, detu) cache stale so the next apply! re-derives from forcing.cmfmc + forcing.dtrain. Called on met-window advance from DrivenSimulation._maybe_advance_window! via the shared invalidate_cmfmc_cache! dispatch.

source

Diffusion

AtmosTransport.Operators.Diffusion Module

julia

Diffusion

Vertical-diffusion operator hierarchy and solver infrastructure.

Public surface:

NoDiffusion — identity no-op; the default when no [diffusion] config section is present.
ImplicitVerticalDiffusion — Backward-Euler implicit diffusion driven by an AbstractTimeVaryingField Kz. Wired into the Strang palindrome via apply_vertical_diffusion_vmr! and installed into TransportModel.diffusion by the runtime recipe when [diffusion] kind = "constant".

Both subtype the global AbstractDiffusion declared in src/Operators/AbstractOperators.jl; concrete operator structs live in operators.jl. The column-level Thomas solve (solve_tridiagonal!) is exposed as the numerical reference. Production kernels name the tridiagonal coefficients (a, b, c) explicitly; the matching transpose lives in src/Adjoints/DiffusionAdjoint.jl.

source

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.

source

AtmosTransport.Operators.Diffusion.DiffusionWorkspace Type

julia

DiffusionWorkspace

Preallocated storage for implicit vertical diffusion.

factors stores the tracer-independent Thomas superdiagonal factors.
layer_thickness stores geometric layer thickness [m] used to convert cell-centered Kz [m² s⁻¹] into interface mass exchange [kg s⁻¹].
references stores one cancellation-reducing column reference per packed cubed-sphere tracer. It is nothing for LL and reduced-Gaussian layouts.

The workspace owns only diffusion data. Advection and convection buffers live in their respective workspace types.

source

AtmosTransport.Operators.Diffusion.DiffusiveSurfaceFluxBoundary Type

julia

DiffusiveSurfaceFluxBoundary()

Inject surface mass immediately before one full implicit vertical diffusion solve, so the new surface mass participates in that solve.

source

AtmosTransport.Operators.Diffusion.ImplicitVerticalDiffusion Type

julia

ImplicitVerticalDiffusion(; kz_field)

Backward-Euler vertical diffusion driven by either a cell-centered Kz field or an exact precomputed TM5 interface dkg 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 with update_field!(op.kz_field, current_time(meteo)); meteo = nothing uses t = 0 for standalone calls.
For Kz fields, reads workspace.layer_thickness 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.
PrecomputedCSDkgField bypasses Kz/geometry reconstruction and does not read layer_thickness; its interface exchange [kg s⁻¹] is already complete.
Uses workspace.factors as Thomas-forward-elimination storage.
Launches a topology-specific mass-flux kernel. Packed layouts factor each atmospheric column once and advance every tracer with those factors.

The spatial operator is linear, but Backward Euler is not a semigroup: V(dt) and V(dt/2) ∘ V(dt/2) differ by O(dt²). The surface-coupling policy therefore selects the timestep composition explicitly.

Fields

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

source

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.

source

AtmosTransport.Operators.Diffusion.SplitSurfaceFluxCoupling Type

julia

SplitSurfaceFluxCoupling()

Apply surface flux at the transport-palindrome center, bracketed by two half-step diffusion solves.

source

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.

source

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.

source

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.

source

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

source

AtmosTransport.Operators.Diffusion.fill_dz_hydrostatic_virtualT! Method

julia

fill_dz_hydrostatic_virtualT!(dz, t_lyr, qv_lyr, ps, ak_ifc, bk_ifc;
                               R = 287.04, gravity = 9.81)

Populate a 3D (Nx, Ny, Nz) dz array using virtual temperature per layer: T_v = T · (1 + 0.61 · qv). Matches the geometry the LocalHoltslagBovilleKzField uses for its column-mid heights, closing the inconsistency between solver dz and Kz-cache dz.

t_lyr, qv_lyr are layer-center 3D fields with shape (Nx, Ny, Nz), typically pulled from the active window's VDIFF payload.

source

AtmosTransport.Operators.Diffusion.fill_dz_hydrostatic_virtualT! Method

julia

fill_dz_hydrostatic_virtualT!(dz_panels::NTuple{6}, t_panels::NTuple{6},
                               qv_panels::NTuple{6}, ps_panels::NTuple{6},
                               ak_ifc, bk_ifc; ...)

Cubed-sphere variant of fill_dz_hydrostatic_virtualT!. Per-panel 3D (Nc, Nc, Nz) dz/t_lyr/qv_lyr arrays + per-panel (Nc, Nc) ps. The VDIFF payload is panel-native, so passing the panel tuples directly is the natural runtime shape.

source

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 (this routine is forward-only)

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

The CS adjoint kernel _vertical_diffusion_cs_single_adjoint_kernel! in src/Adjoints/DiffusionAdjoint.jl builds (a_T, b_T, c_T) inline and runs an inlined Thomas (rather than reusing this routine) so the kernel stays allocation-free.

source

SurfaceFlux

AtmosTransport.Operators.SurfaceFlux Module

julia

SurfaceFlux

Surface emission operators.

Ships the data types and helpers needed to apply per-tracer surface sources to a CellState:

SurfaceFluxSource{RateT} — single-tracer source plus a rate array in model-storage units per second per cell, also available from the top-level AtmosTransport module. Physical inventories in kg species/s are converted to these units by the runtime's basis-aware input builders.
PerTracerFluxMap{S} — NTuple-backed map of SurfaceFluxSources, keyed by tracer_name. Ships with the flux_for(map, :name) lookup helper. Storage-bits-stable on GPU.

The AbstractSurfaceFluxOperator hierarchy (NoSurfaceFlux, SurfaceFluxOperator), the _surface_flux_kernel! KA kernel, and the apply! / apply_surface_flux! entry points live alongside these data types.

Surface layer convention

All kernels here assume k = Nz is the surface. This matches src/Models/DrivenSimulation.jl and the LatLon grid storage layout. A future AbstractLayerOrdering{TopDown, BottomUp} refactor can generalise this.

source

AtmosTransport.Operators.SurfaceFlux.AbstractFluxTemporalScheme Type

julia

AbstractFluxTemporalScheme

How a TimeVaryingSurfaceFluxSource reconstructs a per-step rate from its stored time slices at t = current_time(meteo) (and step dt). Concrete schemes implement

julia

_flux_temporal_weights(scheme, times, t, dt) -> (i0, i1, w0, w1)

returning two bracketing slice indices and blend weights (w0 + w1 == 1) applied as w0·series[i0] + w1·series[i1]. New reconstructions (e.g. mass-flux-integral matching, higher-order in time) are added as new subtypes + one _flux_temporal_weights method — no kernel or operator change.

source

AtmosTransport.Operators.SurfaceFlux.AbstractSurfaceFluxOperator Type

julia

AbstractSurfaceFluxOperator

Top of the surface emission operator hierarchy. Concrete subtypes:

NoSurfaceFlux — identity (default).
SurfaceFluxOperator — wraps a PerTracerFluxMap and applies each tracer's surface flux to k = Nz of the conservative tracer-storage 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 and the reduced-Gaussian face-indexed transport block.

source

AtmosTransport.Operators.SurfaceFlux.AbstractSurfaceFluxSource Type

julia

AbstractSurfaceFluxSource

Supertype for all per-tracer surface emission sources consumed by the surface-flux operator and stored in a PerTracerFluxMap.

Concrete subtypes:

SurfaceFluxSource — a single static per-cell rate array.
TimeVaryingSurfaceFluxSource — a time series of per-cell rate slices plus a sorted times vector; the operator linearly interpolates to the simulation clock at apply! time.

Every subtype carries a tracer_name :: Symbol field used by flux_for to key the map.

source

AtmosTransport.Operators.SurfaceFlux.ConservativeMeanFlux Type

julia

ConservativeMeanFlux()

Window/step integral-conserving: evaluate the linear reconstruction at the step centre t + dt/2. Summed over the substeps of a met window, this midpoint quadrature reproduces the exact trapezoidal time-integral of a piecewise-linear flux — i.e. the emitted mass per window matches the source's own integral, independent of the substep count.

source

AtmosTransport.Operators.SurfaceFlux.LinearInterpFlux Type

julia

LinearInterpFlux()

Linear interpolation between the two bracketing slices, point-evaluated at the step's current_time(meteo). HEMCO-like treatment of instantaneous 3-hourly fields. Default scheme.

source

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.

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 with zero floating-point work and a bit-exact no-op.

source

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

Runtime source construction naturally produces a small tuple; the map keeps that storage and adds 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)    # model-storage units/s per cell
sf6_rate  = fill(1.5e-9, Nx, Ny)    # already converted from physical flux
rn222_rate = fill(3.0e-11, Nx, Ny)

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, <:AbstractSurfaceFluxSource} for some N. Each entry carries a tracer name and its rate data (a static array for SurfaceFluxSource, a time series for TimeVaryingSurfaceFluxSource).

source

AtmosTransport.Operators.SurfaceFlux.PerTracerFluxMap Method

julia

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

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

source

AtmosTransport.Operators.SurfaceFlux.PerTracerFluxMap Method

julia

PerTracerFluxMap(sources::AbstractSurfaceFluxSource...)

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

source

AtmosTransport.Operators.SurfaceFlux.StepwiseFlux Type

julia

StepwiseFlux()

Piecewise-constant in time: emit the value of the slice block containing t (the largest slice time ≤ t), with no blending. This holds each 3-hourly CAMS value constant across its block — the "change the flux in 3h blocks" option.

source

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 conservative tracer-storage 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.

The rate is a model-storage amount per second per cell (already area-integrated); no cell-area multiplier appears in the kernel. File-based physical kg-species rates are converted by the source builder before they reach this operator.

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.

source

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

source

AtmosTransport.Operators.SurfaceFlux.SurfaceFluxSource Type

julia

SurfaceFluxSource{RateT}

A single-tracer surface source: a tracer_name plus a cell_mass_rate array supplying model-storage amount 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 rates are already area-integrated. For dry-VMR tracers, file-based physical fluxes in kg species/s are converted by the source builder to dry-air-equivalent storage units before reaching this struct. The surface flux kernel applies rm_surface += rate × dt without multiplying by cell area.

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

The prognostic tracer is stored per cell, so a per-cell rate × dt is the natural unit. A per-area variant would need to multiply by cell area before entering the operator.

Provenance

The source lives with SurfaceFluxOperator; model and simulation layers only compose and schedule the operator.

Fields

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

source

AtmosTransport.Operators.SurfaceFlux.TimeVaryingSurfaceFluxSource Type

julia

TimeVaryingSurfaceFluxSource{RateT, T}

A single-tracer surface source whose per-cell rate advances through a time series of slices with the simulation clock. Used to carry the CAMS/LMDZ 3-hourly diurnal cycle (instead of collapsing it to a monthly mean) and to match GeosChem.

tracer_name :: Symbol — matches a name in CellState.tracer_names.
cell_mass_rate_series :: RateT — a stack of static rate slices, in the same per-cell model-storage units as SurfaceFluxSource. For cubed-sphere this is an NTuple{6} of (Nc, Nc, ntime) arrays (one stacked panel per face). The trailing dimension is time.
times :: T — slice times in seconds since run start, length ntime, sorted strictly ascending. Stays a host Vector even after Adapt-to-device so the operator can searchsortedlast it on the host before launching the interpolated kernel.
scheme :: S — an AbstractFluxTemporalScheme selecting how the stored slices are reconstructed onto the simulation clock (StepwiseFlux, LinearInterpFlux, or ConservativeMeanFlux). A singleton, carried through Adapt unchanged.

The operator resolves (i0, i1, w0, w1) = _flux_temporal_weights(scheme, times, current_time(meteo), dt) on the host, then blends the two bracketing slices in the kernel. End slices are clamped (constant extrapolation) outside [times[1], times[end]].

source

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

source

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

source

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.

source

AtmosTransport.Operators.SurfaceFlux.flux_temporal_scheme Method

julia

flux_temporal_scheme(name) -> AbstractFluxTemporalScheme

Map a config string to a temporal scheme. "stepwise"/"block" → StepwiseFlux; "linear"/"interp" → LinearInterpFlux; "conservative"/"window_mean"/"integral" → ConservativeMeanFlux.

source

Chemistry

AtmosTransport.Operators.Chemistry Module

julia

Chemistry

Source/sink operators for tracer transformations (decay, photolysis, ...).

Type hierarchy:

julia

AbstractChemistryOperator
├── NoChemistry                     — identity / inert tracers
├── ExponentialDecay{FT, N}         — multi-tracer first-order decay
└── CompositeChemistry              — sequential composition

Interface (OPERATOR_COMPOSITION.md §6):

julia

apply!(state::CellState, meteo, grid, op::AbstractChemistryOperator, dt;
       workspace=nothing)

The operator mutates state.tracers_raw in place and returns state. meteo, grid, and workspace are accepted for interface conformance and may be nothing for operators that do not need them (pure decay).

Multi-tracer decay is fused into a single KernelAbstractions kernel — see chemistry_kernels.jl. Tracers not listed in the operator's tracer_names are left untouched.

source

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

source

AtmosTransport.Operators.Chemistry.ExponentialDecay Type

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

source

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⁻¹

source

AtmosTransport.Operators.Chemistry.NoChemistry Type

julia

NoChemistry()

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

source

Operators API ​

Operators (top-level) ​

Advection ​

Convection ​

Diffusion ​

SurfaceFlux ​

Chemistry ​

Operators API

`Operators` (top-level)

Advection

Convection

Diffusion

SurfaceFlux

Chemistry