Vertical Transport: From Winds to Tracer Columns

This page follows the steps that set vertical tracer transport in AtmosTransport, from met winds to the vertical tracer flux, and compares each step with TM5 and GEOS-Chem High Performance (GCHP). It describes the code as of the CATRINE C90 benchmark (2026-10) and is updated as options are added. Open items are tracked in docs/memos/2026-10-07_v0.5_transport_fidelity_backlog.md.

The vertical index k counts layers from the model top (k = 1) to the surface (k = Nz). Interface k + ½ lies below layer k. Hybrid interfaces are p_{k+½} = A_{k+½} + B_{k+½} p_s, with ΔA_k, ΔB_k the differences across layer k. Above the hybrid region (about 177 hPa for GEOS L72) every layer is pure pressure, ΔB_k = 0.

1. Horizontal mass fluxes

Preprocessing turns met winds into horizontal air-mass fluxes through the cell faces of the cubed sphere, for each layer and substep.

At write time, the flux-fidelity gate independently checks deterministic face samples against the defining integral F = (Δt/g) ∫face [(u eλ + v eφ) ⋅ n] Δp dl. The check runs before column balance, whose deliberate correction is reported separately. Its line-integral evaluator independently rebuilds the 16 midpoint quadrature and four-cell interpolation weights from the meshes instead of applying the bulk sparse operator. Its precision-derived limits are T_abs = 64 eps(FT) max(F_scale, 1 kg) and T_rel = 256 eps(FT); a sample fails only if both its absolute error and its relative error exceed those limits. The check shares its inputs with production, so it cannot detect input errors; its vector evaluator also shares the production face geometry, and native GEOS fluxes are not rederivable (certificate not_applicable). See Binary format for sampling, normalization, and header details.

MERRA-2 (src/Preprocessing/transport_binary/merra2_latlon_regrid.jl).

• Inputs: A3dyn 3-hour mean winds (0.5° × 0.625°); I3 instantaneous surface pressure and humidity.

• u_east and v_north are regridded conservatively to C90 cell means as a Cartesian vector by default ([preprocessing] wind_regrid = "cartesian"); "scalar" retains the historical component-wise path (see Poles below).

• Hourly windows split each 3-hour block; air mass, p_s, humidity and temperature are linear in time.

Two constructions of the face fluxes from cube-centre winds are available ([preprocessing] face_fluxes, MERRA-2 and ERA5 N320 paths); ERA5 N320 also offers line integrals of the source flow ("line_integral", see ERA5 below):

• panel_average (historical, src/Preprocessing/cs_flux_reconstruction.jl, cs_face_fluxes!).

  • Each cell's wind is projected onto its panel's local face-normal directions.

  • The two adjacent cells' projections are averaged.

  • The result is multiplied by the mean moist layer thickness Δp_k, a length, and the time step, and divided by g.

  • The length is the centerline width of the cell on one side of the face (face_lengths = "cell_centerline"). "edge" uses the face's own great-circle length instead.

  • At a panel seam only the cell on the owning panel is used.

• vector (default, cs_vector_face_fluxes!).

  • The two cells' winds are combined as 3-D vectors.

  • At a panel seam the result is further interpolated along the edge to the face midpoint, as FV3 does.

  • It is then projected onto the face's own unit normal and multiplied by the face's great-circle length and the interpolated layer thickness.

Why the seams matter. The cells on both sides of a panel seam are skewed along the edge. They are mirror images across the seam, so both centres are shifted the same way along the edge and the shift does not cancel. The point between the two centres lies 16 % (RMS), up to 25 %, of a face length along the edge from the face midpoint; in the panel interior it is 2e-4. (The two centres are equally far from the face, to 1e-12.) Averaging the two cells therefore takes the wind at the wrong point.

For a solid-body rotation, which has no divergence, the spurious divergence relative to the face flux at C90 is as below (/temp1/cfranken/scratch/overnight_2026_10_07/solid_body_divergence.jl; test/core/test_cs_face_fluxes.jl checks the same at C48):

	interior cells	seam cells
panel_average	1e-4	4e-3
vector	2e-7	5e-6

In pure-pressure layers nothing corrects this divergence (section 2), so it becomes vertical motion. Against MERRA-2's own vertical velocity (2022-07-15, C90, p < 150 hPa), the panel_average error at seam cells is larger than the true signal: RMSE 8.8e-3 Pa/s against an RMS of 7.2e-3. With vector it is 1.3e-3, as in the interior. In the extratropics at 30–200 hPa the RMSE falls from 8.4e-3 to 4.6e-3 Pa/s.

In a four-month tracer run against GCHP (December 2021 – March 2022, ours started from GCHP's first output, hybrid_mass and the FV3 vertical profile), vector changes:

• the CO₂ growth of the bias above 100 hPa from +0.050 to −0.035 ppm, and the SF₆ growth from +0.006 to −0.004 ppt;

• the slope against GCHP above 100 hPa from 0.979 to 0.998;

• the near-surface Rn-222 RMSE by −5 % and the CO₂ RMSE from 0.331 to 0.310 ppm.

With in addition wind_regrid = "cartesian" and face_interpolation = "fv3" (GCHP's construction, now the default):

• the stratospheric growth is the same (−0.036 ppm);

• the near-surface RMSE falls further, by 15 % (CO₂), 16 % (fossil CO₂) and 24 % (Rn-222), and the column CO₂ RMSE from 0.101 to 0.085 ppm.

"cubic" without the filter gives mixed near-surface changes (worse for CO₂, fossil CO₂ and SF₆, better for Rn-222).

The ERA5 N320 and spectral-to-CS paths (era5_n320_regrid.jl and cubed_sphere_spectral.jl) take the same options and defaults. The separate lat-lon-binary regrid path (cubed_sphere_regrid.jl) retains the historical construction because its input binary does not preserve endpoint source winds. The reduced Gaussian source cells have great-circle edges along their latitude bounds, so they leave up to 0.3 % of the polar C90 cells uncovered. The N320 path therefore divides every intensive field by the regridded constant 1 (_regrid_intensive!). Before 2026-10-08 a constant field regridded to 0.9967 poleward of 89°: surface pressure was about 3 hPa low and temperature about 0.8 K cold in the 8 polar cells.

Poles. The east and north unit vectors turn with longitude. Regridding u_east and v_north as two scalars therefore averages components that point in different directions. Each of the four C90 cells that meet at a pole spans 90° of longitude. MERRA-2's pole row holds one wind vector, written in the local basis of each longitude. Its Cartesian components vary by only 6–8% along the row, while u itself changes sign.

Effect of scalar regridding:

• On an analytic flow across the pole it puts errors of about 5% in the wind poleward of 88°, and 0.5% at 85–88°.

• After vector face fluxes, the resulting spurious divergence at |lat| > 88° is 2.8 times the global RMS divergence, and 0.28 times it at 85–88° (/temp1/cfranken/scratch/overnight_2026_10_07/polar_wind_regrid.jl).

• Against MERRA-2's vertical velocity (2022-07-15, p < 150 hPa), the ω error poleward of 88° is 7.4e-3 Pa/s, against a true RMS of 2.1e-3 Pa/s. At 85–88° it is 1.0e-3 against 2.1e-3.

wind_regrid = "cartesian" (CartesianWindRegrid) does what MAPL does for GCHP's UA;VA:

• writes the source winds as Cartesian components;

• regrids those three components conservatively;

• projects the result back onto east and north at each cube-cell centre.

GCHP's own C90 winds (Met_U/Met_V in the CATRINE output) match it to within 5e-4 m/s at every latitude and level. They differ from the scalar regrid by 0.5–1.5 m/s poleward of 88° (/temp1/cfranken/scratch/overnight_2026_10_07/gchp_wind_compare.jl).

Interpolation to the faces. face_interpolation (vector only) chooses how the cell winds reach the face midpoint:

• "linear": the two adjacent cells.

• "cubic": FV3's fourth-order stencil (−1, 9, 9, −1)/16 across interior faces that have two cells on each side in the same panel row. The layer thickness stays linear. Against exact line integrals of a smooth analytic wind (wavelength about 30 cells), interior face fluxes are three times more accurate (face_flux_accuracy.jl). The divergence gains less, 1.6 times at that wavelength, because the faces next to the seams stay second order (CSVectorFaceGeometry); for very smooth winds it is slightly worse.

• "fv3" (default): "cubic" plus the filter that GCHP's restaggering applies along each face (next paragraph).

flux_thickness = "dry_mass" uses the dry layer thickness that the fluxes transport, instead of total-air Δp:

$$\mathrm{\Delta }{p}_{\mathrm{dry},k}=(1-q_k)|\mathrm{\Delta }{A}_k+\mathrm{\Delta }{B}_k{p}_s|=\frac{g{m}_{\mathrm{dry},k}}{A_{\mathrm{cell}}}.$$

For target-grid face methods this thickness is evaluated in each cubed-sphere cell before interpolation to faces. ERA5's line-integral method evaluates the equivalent (1-q) Δp on N320 before integrating along a cube face. ERA5 N320 defaults to dry_mass, matching its dry-basis layer masses; MERRA-2 retains its historical moist default for compatibility.

ERA5 (.../era5_n320_regrid.jl and .../cubed_sphere_spectral.jl). Instantaneous hourly winds are supplied on the N320 reduced Gaussian grid or synthesized onto the regular lat–lon staging grid. Surface pressure comes from ARCO or the spectral LNSP field.

• With panel_average or vector the winds are regridded to the cube as above.

• face_fluxes = "line_integral" (LineIntegralFaceFluxes) integrates the source-grid flow along each cube face instead:

  $$F_f=\frac{\mathrm{\Delta }t}{g}{\int }_f(\mathbf{V}\cdot {\mathbf{N}}_f)\mathrm{\Delta }pdl.$$

  A cube face is a great-circle arc, so its unit normal N_f is the same vector along it. The Cartesian components of V Δp are interpolated bilinearly from the N320 or staging-grid cells to 16 midpoints of the arc. By the divergence theorem the net flux out of a cube cell is then the N320 flow's divergence integrated over the cell; the cell-centre methods smooth it at the cube grid scale. TM5 integrates the spectral winds along its cell edges in the same spirit.

• flux_time_sampling = "window_mean" (default) uses (F(h) + F(h + 1)) / 2, the trapezoidal rule for the window's mean flux; the Poisson balance then closes it against m(h + 1) − m(h) as before. "window_start" holds F(h) over the window. Without a next day on disk the last window uses F(h) and records the boundary fallback.

Against ERA5's own vertical velocity (ARCO-ERA5 ar/model-level-1h-0p25deg.zarr-v1, 2022-01-15, C90, hybrid_b), in the pure-pressure layers above 71 hPa, where ω = g cm / A:

face fluxes	RMSE / RMS of ERA5 ω	slope per C90 cell	slope on 3 × 3 blocks
vector, fv3 filter, L137	0.37–0.39	0.73–0.75	0.92–0.94
line_integral, L117	0.03 (10–71 hPa), 0.15 (< 10 hPa)	0.95–1.00	0.94–0.99

Both are compared with ω at the window start, with window_start sampling; with window_mean the line-integral fluxes match the mean of ω at the start and the end of the window equally well. Above 10 hPa the L117 layers merge several native levels, and the mean of the interface fluxes differs from the layer mean of ω. In hybrid layers above 200 hPa the slope is 1.01 (0.74 with vector). The line integrals thus resolve ERA5's vertical motion at the cube grid scale. Its larger grid-scale values raise the adaptive substep count: 11–20 per window (mean 14) with line_integral and window_mean on L117, against 10–13 with vector on L137 and 9–10 with vector on L66. The largest column mismatch before the Poisson balance falls from 2.6e11 to 1.2e11 kg. Scripts and output: /temp1/cfranken/scratch/arco_ml_check/.

GCHP takes the same MERRA-2 A3dyn winds. MAPL regrids them conservatively to the cube as a vector (three Cartesian components). GCHPctmEnv then restaggers them with A2D2C (GEOS_FV3_Utilities.F90, d2a2c_vect in sw_core.F90):

1. A → D. The along-edge wind at each face is the average of the two cells' Cartesian winds. At cube edges it is interpolated along the edge (edge_vect_*).

2. D → A, fourth order (dord4). In the panel interior (cells at least four from an edge) the step is (−1, 9, 9, −1)/16; nearer the edges it is the two-point mean. Together with step 1 this returns the A-grid wind filtered along the face direction: (−1, 8, 18, 8, −1)/32 in the interior and (1, 2, 1)/4 near the edges. The filter removes the 2Δ wave along the face and keeps 62% of the 4Δ wave.

3. A → C. Fourth order (−1, 9, 9, −1)/16 across the face in the interior, one-sided next to the edges, and edge_interpolate4 on the edges.

fv_computeMassFluxes then forms contravariant winds and the swept area u dt dy sin_sg. It advects the start-of-step dry Δp = ΔA + ΔB p_s,dry with FV3's PPM over nsplt subcycles (fv_tp_2d), so the face Δp is an upwind-swept value, not the two-cell mean. In pure-pressure layers Δp is uniform at the start, so the first subcycle agrees with the two-cell mean; later subcycles see the Δp updated by the divergence.

face_interpolation = "fv3" reproduces the interior stencils of steps 1–3. Four differences from GCHP remain:

• FV3 interpolates covariant components and converts them to contravariant ones with the grid angles; we project the vector on the face normal.

• FV3 picks the D → A order per contributing cell; we pick the filter per face from the face's position. They differ in a band three to five cells from the panel edges.

• The first and last row of a panel are not filtered here; FV3 filters them using halo cells from the neighbouring panel.

• Next to the edges FV3 uses its one-sided and edge_interpolate4 interpolations; we use the two-cell value with the along-edge correction.

The winds are held over each 3-hour block. GCHP reads the A3dyn record at t + 90 min, so record r advects hours 3(r−1) to 3r, as here. The CATRINE output at 00z holds the previous day's last record: it is the step that ends at 00z. PS1/PS2 are interpolated in time to each 600 s step (GCHP run-directory setting IMPORT_MASS_FLUX_FROM_EXTDATA = .false. for MERRA-2).

TM5 integrates U·Δp along each cell edge from the spectral fields (tmm.F90), with at least two points per spectral wavelength.

2. Column mass closure (pressure fixer)

The reconstructed fluxes never close the column mass budget exactly. The mismatch between flux convergence and the met surface-pressure tendency is about 0.1–0.2 Pa/s RMS before correction. AtmosTransport closes each column:

1. Solve a Poisson problem on the cube for a correction to the column-integrated face fluxes, so that column convergence equals the change of column air mass over the window (balance_cs_column_mass_fluxes!, src/Preprocessing/cs_poisson_balance.jl).

2. Spread each face's column correction δF over the levels with weights w_k: δF_k = w_k δF / Σ_j w_j (spread_column_correction!).

The lat-lon (balance_column_mass_fluxes!) and reduced-Gaussian (balance_compressed_column_fluxes!) spectral paths close their columns the same way, with the same weights. The explicit-tendency cm closure (recompute_cm_from_dm_target!, diagnose_cs_cm! on the cube) folds a column's remaining residual R into the levels with the same weights, C_{k+1} ← C_{k+1} − R Σ_{j≤k} w_j / Σ_j w_j. No path balances each level to its own mass change: with the column closed, the per-level imbalance that remains is the vertical mass flux of section 3, and a per-level balance would remove it. TM5's per-level balance (BalanceMassFluxes_m) keeps vertical motion only because its target includes the IFS vertical flux; our preprocessing has no independent vertical flux. The reduced-Gaussian path balanced per level until issue #11, which left its vertical mass flux at about 1e-4 of the lat-lon one on the same ERA5 day.

The weights are selectable (column_balance_weights, on every path that applies a column balance: lat-lon, reduced Gaussian and cubed sphere; Data sources lists where each path sets it):

Option	Weights w_k	Consequence
mass	layer air mass of the two adjacent cells	The correction reaches every layer, including the stratosphere
hybrid_b	ΔB_k	Only hybrid layers are corrected; the semantics of TM5 and GCHP
hybrid_mass (default)	layer air mass where ΔB_k ≥ 10⁻⁶, else 0	Like mass in the troposphere, no correction in pure-pressure layers

hybrid_mass counts a layer as hybrid when ΔB_k ≥ 10⁻⁶ (HYBRID_LAYER_MIN_DB): its thickness Δp_k = ΔA_k + ΔB_k p_s then changes by at least 10⁻⁶ × 10⁵ Pa = 0.1 Pa when p_s changes by 1000 hPa. GEOS/MERRA-2 L72 layer 41 (150–177 hPa) has ΔB = 8.2 × 10⁻⁹ (8.2 × 10⁻⁴ Pa per 1000 hPa), a rounding residue of the B table (3 × 10⁻⁷ of its 2654 Pa thickness); with a plain ΔB > 0 mask it took about 3 % of every column correction at full mass weight although its thickness does not follow p_s. The smallest genuine hybrid increment of the shipped tables is ERA5 L137 layer 55 (73–77 hPa, the topmost layer with ΔB > 0: ΔB = 7 × 10⁻⁶, 0.7 Pa per 1000 hPa, 0.17 % of its 416 Pa thickness at p_s = 1000 hPa); the ERA5 level presets keep it at ΔB = 7 × 10⁻⁶ (ml137_full, ml137_cfl94, ml137_cfl85) or merge it into a layer with 2.4 × 10⁻⁵ (ml137_66L) or 5.9 × 10⁻⁵ (ml137_tropo34, ml137_tropo34a), all above the threshold. The next GEOS layer has ΔB = 7 × 10⁻³ (696 Pa per 1000 hPa). On the GEOS-IT C45 golden day the threshold changes the layer-41 fluxes by 1.8–1.9 % RMS and cm at 177–288 hPa by 12 % to 5 % RMS; on MERRA-2 C90 by 0.25 % and 3.5 % to 1.4 %, and the MERRA-2 OMEGA NRMSE improves from 0.4725 to 0.4712.

Size of the default change. hybrid_mass against mass on the golden days (2021-12-01 GEOS-IT and ERA5, 2022-01-15 MERRA-2; the cubed-sphere numbers were measured with the ΔB > 0 mask, before the threshold above):

• GEOS-IT native C180 → C45: stratospheric RMS cm falls by 25–43 % between 5 and 130 hPa (0.57× at 48 hPa). Horizontal fluxes on faces above 10 % of the global maximum change by up to 17.8 % (99th percentile 2 %); flux-weighted, by 1–2.7 % in the pure-pressure layers 1–40 and 0.25–0.7 % in the hybrid layers. Tropospheric cm changes by less than 0.3 %. With mass weights the column residual is pushed into the pure-pressure layers as spurious stratospheric vertical flux (a young-air bias); hybrid_mass keeps their native divergence, as GCHP's remap in fixed-pressure layers and TM5's ΔB weighting do.

• ERA5 spectral C24: stratospheric RMS cm 0.65–0.72×.

• ERA5 lat-lon 72×37: stratospheric RMS cm 0.60–0.69× at 10–100 hPa and 0.72–0.93× at 100–200 hPa; reduced Gaussian O24: 0.70–0.83× and 0.83–0.96×.

• Air masses, ps and dm are unchanged in every case.

The vertical mass flux (section 3) inherits the correction as Σ_{j ≤ k} (w_j/Σw − ΔB_j) · δ, a mode that is coherent in the vertical. With mass weights it is proportional to p/p_s − B, about 0.1 of the column mismatch at 100 hPa and 0.2 near the extratropical tropopause.

Measured effect. These tests compare the vertical velocity diagnosed from cm with MERRA-2's own (A3dyn OMEGA, C90, 3-hour means, on 2022-01-15 and 2022-07-15).

• At 30–100 hPa poleward of 30°, the RMS of the mass-weighted vertical velocity is 1.0–3.5 times that of OMEGA. The excess is largest at high summer latitudes, where the resolved motion is weakest. The correlation is 0.28–0.93. With hybrid_mass the RMS ratio is 0.97–1.10 and the correlation 0.87–0.99.

• Below 200 hPa, hybrid_mass and mass agree with OMEGA equally well (correlations within 0.01).

• hybrid_b matches hybrid_mass in the stratosphere. In the troposphere it is close but slightly worse, mostly in the southern storm track (60–30°S). The RMS difference from OMEGA, in 1e-3 Pa/s, is:

  • at 500–850 hPa: 69.9 against 65.8 (January), 76.5 against 69.4 (July);

  • at 200–500 hPa: 59.4 against 52.3 (January), 65.7 against 56.2 (July).

Per-band values: /temp1/cfranken/scratch/omega_3way_2022{0115,0715}_dB.txt (memo 2026-10-07_v0.5_transport_fidelity_backlog.md, item 1). A first hybrid_b build had a constructor bug: its weights were the B value at each layer's top interface, not ΔB. Its results are superseded by these.

In a four-month MERRA-2 run (December 2021 – March 2022, PPM, compared with GCHP), hybrid_mass cut the growth of the bias above 100 hPa from +0.54 to +0.19 ppm CO₂ and from +0.066 to +0.023 ppt SF₆. The slope against GCHP above 100 hPa rose from 0.89 to 0.97. Near-surface RMSE (below 910 hPa) changed little:

• CO₂ fell from 0.383 to 0.370 ppm and SF₆ by 5 %;

• fossil CO₂ and Rn-222 rose by 0.3–0.4 %.

Column RMSE was unchanged. hybrid_b gives the same stratospheric improvement (+0.19 ppm CO₂) with near-surface RMSE between the two: CO₂ 0.377 ppm.

These growth values include the decay of an initial difference.

• GCHP's initial lower stratosphere lies 0.5 ppm CO₂ and 0.06 ppt SF₆ above the protocol initial state at 50–100 hPa; ours matches the protocol file. GCHP's state is reproduced, to 0.05 ppm and 0.003 ppt, by placing each protocol layer value at the layer's top interface instead of its middle.

• Started from GCHP's own first output instead, hybrid_mass with the FV3 vertical profile grows by +0.05 ppm CO₂ above 100 hPa in four months, not +0.14 (memo 2026-10-07_v0.5_transport_fidelity_backlog.md, section B1).

• Differences between options run from the same initial state are not affected.

GCHP has no column closure. FV3 transports tracers horizontally with the uncorrected fluxes. It then remaps the Lagrangian layers onto A + B p_s,adv, where p_s,adv is the surface pressure the flux convergence implies (fv_tracer2d.F90, offline_tracer_advection). The mismatch with the met surface pressure stays in p_s,adv. The pressure is reset from the met fields at the next step. Each tracer's global mass is then rescaled, with the factor applied only from the first layer with B > 0 downward (fv_tracer2d.F90, about line 977). The implied vertical motion distributes the column convergence by ΔB.

TM5 holds its vertical flux at the IFS value and Poisson-corrects the horizontal fluxes level by level. The column tendency it closes is ΔB_k (p_s2 − p_s1), so the mismatch never reaches pure-pressure layers.

3. Vertical mass flux

The vertical air-mass flux through interface k + ½ follows from continuity for each substep (diagnose_cs_cm!):

cm_{k+½} = cm_{k−½} + (∇·F)_k − Δm_k, with cm_½ = 0 at the top.

Δm_k is the window's change of layer air mass per substep. A residual left at the surface interface is spread with the same weights as in section 2. The binary stores cm for each window; it is constant within the window. GCHP has no stored vertical flux: its vertical motion is implicit in the remap. TM5 recomputes cm at every dynamics step from the time-interpolated fluxes.

When source OMEGA is available, write-time certification additionally compares it with ω_cm = g cm/(A_cell Δt_amount) after matching dry (1-q) basis, time sampling, horizontal grid, and vertical levels. The header reports global and pressure-band errors against a named per-source baseline, which passes only with NRMSE < 1 (the score of cm ≡ 0), correlation > 0.5 and 0.5 < slope < 1.5 (MERRA-2 C90: 0.47, 0.89, 0.87; the slope bounds reject a cm that is uniformly a factor 2 or more too small or too large). This is a physical-fidelity diagnostic distinct from replay closure: replay asks whether the stored arrays conserve mass together, while OMEGA asks whether their vertical motion agrees with an independent meteorological field.

4. Vertical advection

Splitting. Each substep applies X Y Z [diffusion, emissions] Z Y X (src/Operators/Advection/CubedSphereStrang.jl). The C90 MERRA-2 runs of December 2021 – March 2022 used 6–13 substeps per hourly window (mostly 8–11), that is 12–26 vertical sweeps per hour. Convection and chemistry run once per window. Air mass is updated in every sweep with the same fluxes that move the tracer, so tracer mass is conserved exactly.

Flux form. In each vertical sweep rm_k ← rm_k + f_{k−½} − f_{k+½}, with f = α (rm_donor + (1 − α) s_donor) for downward flow (and the mirror form for upward flow), where α = cm / m_donor and s is a tracer-mass slope.

scheme = "ppm" (production). The slope comes from the PPM edge value on the outflow side of the donor cell:

• edge values use the uniform-index fourth-order formula q_{k+½} = 7/12 (q_k + q_{k+1}) − 1/12 (q_{k−1} + q_{k+2}), without layer-thickness weighting;

• the edges are limited with the Colella–Woodward (1984) monotonicity conditions (flattening at extrema, overshoot correction);

• the limited edge sets a linear slope toward the outflow face, which is bounded by the cell content (_limited_moment), and the flux uses the Russell–Lerner form above (_slopes_face_flux, src/Operators/Advection/reconstruction.jl).

The scheme is therefore second order, with a PPM-informed slope; the parabola is never integrated. The top two and bottom two layers fall back to upwind. In 3-D tests it behaves almost like scheme = "slopes" (minmod slopes).

scheme = "ppm", vertical = "fv3_kord8" (cubed sphere) uses the profile GCHP uses for tracers, FV3's scalar_profile with kord = 8 and iv = 0 (fv_mapz.F90), in the same flux form (src/Operators/Advection/vertical_fv3_profile.jl). This is also the default vertical remap for configured scheme = "linrood" runs, following the runtime rule that the default should be the closest GCHP match. vertical = "upwind" remains available explicitly. Both have an adjoint (see the adjoint of the FV3 vertical sweep). Within layer k, with s ∈ [0, 1] running from the top to the bottom edge,

$$q(s)=q_L+s[(q_R-q_L)+q_6(1-s)],q_6=3(2{\overline{q}}_k-q_L-q_R).$$

1. Edges. The edge values solve FV3's compact tridiagonal system. It is fourth-order accurate on a grid whose spacing is the layer air mass, so thin and thick layers are weighted correctly.

2. Large-scale constraints. Each interior edge lies between its two neighbouring layer means, except at a local extremum of the means. There it may overshoot on the extremum side only, and at a minimum it stays non-negative.

3. Layer limiters.

• Layers 3 to Nz − 2: Huynh's second constraint bounds each edge using the neighbouring slopes, then a positive-definite limiter removes negative interior minima.

• Layers 2 and Nz − 1: the standard PPM monotonicity limiter; at a local extremum of the means the layer is flattened.

• Top and bottom layers: the outer edge is made non-negative, then the layer is monotone.

4. Flux. The tracer flux through an interface is the air-mass flux F times the mean of the donor parabola over the swept fraction α = |F| / m_donor ≤ 1:

$${\overline{q}}_{\downarrow }=q_R-\frac{\alpha }{2}[(q_R-q_L)-(1-\frac{2\alpha }{3})q_6],{\overline{q}}_{\uparrow }=q_L+\frac{\alpha }{2}[(q_R-q_L)+(1-\frac{2\alpha }{3})q_6],$$

for downward flow (bottom of the layer above) and upward flow (top of the layer below).

In cumulative air-mass coordinates, this sweep is FV3's conservative remap of the layers (mapn_tracer, which uses the positive-definite profile) onto interfaces shifted by the swept mass, M_e → M_e − F_e. In exact arithmetic the two are identical, provided that:

• no interface flux exceeds its donor layer's air mass;

• the fluxes through the model top and the surface are zero.

One sweep is one remap with the same profile. GCHP remaps once per 600 s after its horizontal step and then rescales tracer mass globally. fillz is off by default (fill = .false. in fv_arrays.F90). Here the remaps follow the binary's cm, twice per substep (per palindrome subcycle).

Each interface flux is computed once per column and enters both adjacent layers with opposite signs, so the fluxes cancel exactly in the column sum. Only the rounding of the cell updates remains, in Float32 as in Float64.

The profile is non-negative everywhere when the layer means are. A layer's outflow is then the integral of its own non-negative parabola over the swept parts. The sweep therefore stays non-negative as long as the swept fractions of a layer's two faces add up to at most one, α_top + α_bottom ≤ 1. The runtime's subcycling budgets the total outflow of all six sweeps of a substep against the starting layer mass (cfl_limit = 0.95). That budget is a static proxy, not a proof for the evolving mass, but it keeps the sum below one in practice.

Signed tracers. The profile above is FV3's positive-definite one (iv = 0), which GCHP uses for all tracers. It flattens layers whose mean is not positive and clamps edges at zero, so tracers that become negative, such as flux anomalies, fall back to first order where they are negative. vertical = "fv3_kord8_signed" selects FV3's profile for signed fields (iv = 1). It omits the three non-negativity steps and is symmetric under q → −q, except where neighbouring layer means are exactly equal: FV3's large-scale constraint resolves those ties with its local-minimum branch whatever the sign.

Differences from GCHP:

• GCHP remaps Lagrangian layers once per 600 s step; here the profile drives each flux-form vertical sweep, with interface fluxes from the binary's cm.

• GCHP rescales each tracer's global mass in the hybrid layers after every step; this model conserves tracer mass without rescaling.

Tests (test/core/test_fv3_vertical_profile.jl):

• The edge values and limited parabolas equal, bit for bit, an independent transcription of fv_mapz.F90 for both profiles; the positive-definite parabolas are non-negative inside every layer.

• Beyond a Courant number of one the flux is capped at the donor's content.

• The fluxes equal the swept integrals of that reference. The whole update equals an independent remap onto the shifted interfaces, to 1e-11.

• Uniform fields are preserved, column mass telescopes and results stay non-negative, in both precisions.

• A six-panel sweep equals the column kernel column by column, the single-tracer and packed sweeps agree, and TransportModel builds the column scratch and runs the sweep.

Checks on an L40S (C90, L72, four tracers, 200–400 random sweeps at Courant ≤ 0.3; scripts and output in /temp1/cfranken/scratch/fv3_vertical/):

• CPU and GPU air mass are bit-identical.

• Relative tracer-mass changes are below 3e-15 in Float64 and 5e-9 in Float32. Under the same Float32 test the default PPM changes the uniform background tracer by up to 1.4e-8.

• One-month C90 MERRA-2 run (December 2021, production Courant numbers, hybrid_mass binaries):

  • Float32 − Float64 global tracer totals: −3e-8 (CO₂), +6e-8 (SF₆) and −2e-10 (fossil CO₂).

  • The default scheme, from the 25-month twin on mass binaries (so the closure differs as well), gives −4e-8, +5e-8 and +8e-9 for the same month.

  • Column-mean differences match (RMS 0.009 ppm CO₂ for both; /temp1/cfranken/catrine_protocol/compare_f32_f64_dec2021/).

• A six-panel sweep costs 1.2 ms in Float32 against 0.47 ms for the default. In a four-month C90 run the wall time per window is unchanged (0.58–0.60 s).

scheme = "linrood" uses FV3-style cross-term advection horizontally and, by default, first-order upwind vertically. That is far too diffusive: in a four-month MERRA-2 test the stratospheric CO₂ bias against GCHP grew six times faster than with ppm. With vertical = "fv3_kord8" its vertical sweeps use the FV3 profile above.

GCHP has no Eulerian vertical advection. After the horizontal step it remaps the Lagrangian layers onto the hybrid levels with PPM (kord_tr = 8: thickness-weighted edges, Huynh constraint, positivity), once per 600 s step. TM5 carries prognostic tracer-mass slopes that are advected with the tracer, with a positivity-only limiter and xyz–zyx splitting.

Numerical diffusion. In 1-D column tests on the L72 grid (/temp1/cfranken/scratch/tm5_advection/results_matrix.txt), a 5-km vertical wave is advected at ω = 0.0073 Pa/s. That is close to the mean MERRA-2 |ω|; the median per-sweep Courant number is 8e-4. The effective diffusivities at 100–30 hPa are 0.57 m²/s for ppm, 0.30 for TM5's prognostic slopes, 0.08 for FV3 kord 8 (Lagrangian remap and our flux form alike) and 2.7 for upwind. For a smooth monotone profile, ppm and FV3 are near zero, while both slopes schemes and upwind are not.

In 1-D translation tests (a Gaussian bump moved 30 layers; /temp1/cfranken/scratch/fv3_vertical/translation_probe.txt):

• The error of the FV3 profile falls 11–23× when the bump width doubles (third order or better). It falls 4× for ppm and slopes (second order) and 2× for upwind.

• The FV3 profile is 4–900 times more accurate than ppm for widths of 3–24 layers and Courant numbers of 0.01–0.3.

In four-month MERRA-2 runs against GCHP (December 2021 – March 2022), the FV3 profile reduced the growth of the CO₂ bias above 100 hPa by about 0.05 ppm, with either column-balance weighting:

• with mass weights, from +0.54 to +0.49 ppm;

• with hybrid_mass weights, from +0.19 to +0.14 ppm.

The FV3 profile also lowered near-surface RMSE for every tracer, by 1–3 %. The column closure (section 2) is the larger of the two effects.

Adjoint of the FV3 vertical sweep

The cubed-sphere adjoint (cs_surface_emission_footprint, the 4D-Var driver) transposes Lin–Rood's vertical sweeps with the FV3 profile (vertical = "fv3_kord8" or "fv3_kord8_signed") as well as the upwind ones. Meteorology is fixed, so the adjoint carries only the tracer mass.

The sweep as a matrix. In one column, with m the air mass entering the sweep and F_e the interface fluxes, the sweep is

$$\overline{q}=Srm,e=T\overline{q},(q_L,q_R,q_6{)}_k=P_k(\beta )x_k,f_e=F_e^{d}c({\alpha }_e)\cdot (q_L,q_R,q_6{)}_{\mathrm{donor}(e)},r{m}_k^{\prime }=r{m}_k+f_k-f_{k+1},$$

with S = diag(1/m) (here and below m is floored at the advection mass floor, as in the forward), T the compact edge interpolation (a Thomas solve whose coefficients depend on m only), x_k = (e_k, e_{k+1}, q̄_{k−2}, …, q̄_{k+2}) the inputs of layer k's parabola, and F^d, α = |F^d| / m_donor the flux capped at the donor's mass and its swept fraction. Integrating the parabola gives the weights

$$c_{\downarrow }(\alpha )=(\frac{\alpha }{2},1-\frac{\alpha }{2},\frac{\alpha }{2}(1-\frac{2\alpha }{3})),c_{\uparrow }(\alpha )=(1-\frac{\alpha }{2},\frac{\alpha }{2},\frac{\alpha }{2}(1-\frac{2\alpha }{3}))$$

acting on the donor's (q_L, q_R, q_6): a downward flux sweeps the bottom of the layer above the interface, an upward flux the top of the layer below. Every limiter step of section 4 replaces an edge or a parabola by one of a few linear forms of x_k (an edge clipped to a neighbouring mean, 3 q̄ − 2 q_R, a flattened layer (q̄, q̄, 0), zero); comparisons choose which. For fixed choices β the 3 × 7 matrix P_k(β) is constant and the sweep is linear, rm' = W(β) rm.

Reverse pass. The adjoint is the transpose at the forward's branches,

$${\lambda }_{rm}={\lambda }^{\prime }+S(T^{\mathsf{T}}{\lambda }_e+{\lambda }_{\overline{q}}),({\lambda }_e,{\lambda }_{\overline{q}})=\sum _k{P}_k(\beta {)}^{\mathsf{T}}\sum _{e:\mathrm{donor}(e)=k}{F}_e^{d}c({\alpha }_e)({\lambda }_e^{\prime }-{\lambda }_{e-1}^{\prime }),$$

computed column by column (_cs_zsweep_fv3_adjoint_kernel!, src/Operators/Advection/vertical_fv3_profile_adjoint.jl): the interface adjoints λ'_e − λ'_{e−1}, through each donor's parabola to its edge and mean adjoints, then the transposed Thomas solve (back substitution first, then the elimination rows from the surface up), then S.

What is recorded. Each vertical-sweep record holds the air and tracer mass entering the sweep, with cm and the flux scale. The reverse pass replays the forward profile code (_fv3_edge, _fv3_layer_profile, _fv3_interface_flux) on numbers that carry their partial derivatives with respect to x_k. Their value arithmetic is the forward's floating-point arithmetic, so every comparison, and therefore β, is the forward's; the partials give P_k(β) and c(α) without a second copy of the limiter logic. The record adds one full-panel tracer tuple per vertical sweep to the Lin–Rood tape: twelve instead of ten full-panel tuples per step, plus six face tuples. At C90 L72 in Float64 that is 552 MB per step instead of 488 MB (L137: 1050 MB instead of 929 MB). Storing one byte of branch flags per cell instead would save about 11 % of the tape, at the cost of a second implementation of the limiter decisions.

Composition. Per step the Lin–Rood tape is H Z [diffusion] emissions [diffusion] Z H [convection], and the reverse pass applies the transposes in the opposite order. Stride and Revolve checkpointing store (rm, m) and record each window again from its checkpoint, so they reproduce the same branches; their footprints equal those of the full tape to rounding (rtol = 1e-10).

Frozen-limiter caveats. The adjoint is the gradient of the piecewise-linear forward map within the branch region of the recorded state. On a branch boundary (a tie between two candidate values) it is the derivative of the branch the forward took. The PPM monotonicity and positivity limiters switch discontinuously, so the cost is only piecewise smooth. With a zero background, the positive-definite profile sits on its q̄ ≤ 0 branch (a flat layer) in every layer the tracer has not reached. There the adjoint is that of the flattened profile, not of the parabola that a small positive tracer would get. Finite-difference checks therefore need a smooth, nonzero base field, as for the horizontal Lin–Rood adjoint.

Tests (test/core/test_cs_linrood_fv3_adjoint.jl):

• The replayed parabolas equal the forward's bit for bit, on columns that engage every limiter, for both profiles.

• In 256-bit arithmetic, ⟨λ', W δ⟩ from central differences (h = 1e-40) equals ⟨Wᵀ λ', δ⟩ to 1e-30 on such columns, with fluxes in both directions and fluxes capped beyond Courant number 1.

• Lin–Rood footprints match central differences to 1e-9 (observed ≤ 1e-10) for ORD 5 and 7, upward and downward fluxes, surface-layer and column-mean objectives, and both profiles (the signed one on a field that crosses zero); a linearity witness confirms that no limiter switches branch over the stencil.

• Full, Stride and Revolve checkpointing give the same footprints to rounding, and the tape's forward replay equals production transport bit for bit.

The split-sweep scheme = "ppm" with vertical = "fv3_kord8" has no adjoint: its x, y and z records share one scheme type, so the reverse pass cannot yet tell its vertical records apart.

5. Sub-grid vertical transport

Convection and boundary-layer mixing add vertical transport that the met winds do not resolve. See the convection (CMFMC with the DQRCU cloud base; TM5 matrix convection) and diffusion (GEOS-Chem non-local VDIFF; TM5 dkg) documentation. In the benchmark, ERA5 runs use the TM5 schemes and MERRA-2 runs the GEOS-Chem schemes.

6. Summary

Step	TM5	GCHP 14.7 (MERRA-2)	AtmosTransport
Horizontal fluxes	Spectral edge integrals	FV3 from conservatively regridded A3dyn winds	Two-cell means of conservatively regridded winds × layer thickness
Column closure	Per-level Poisson, vertical flux fixed, ΔB	None; mismatch in advected p_s, reset each step	Column Poisson; weights hybrid_mass (default), hybrid_b, mass
Vertical motion	IFS vertical velocity	Implied by Lagrangian remap (ΔB)	cm from continuity, constant per window
Vertical scheme	Prognostic slopes, positivity limiter	Lagrangian remap, PPM kord 8	Flux form; PPM-informed slope with CW84 limiter (ppm), or FV3 kord 8 parabola (vertical = "fv3_kord8")
Time step	Global ndyn, CFL-limited	600 s	10–13 substeps per hourly window

GEOS-Chem is a reference, not the physical truth. The choices above are judged against the met model's own diagnostics and against observations, as well as against GCHP.

References

• Colella, P., and Woodward, P. R. (1984). The piecewise parabolic method (PPM) for gas-dynamical simulations. J. Comput. Phys. 54, 174–201.

• Eastham, S. D., et al. (2018). GEOS-Chem High Performance (GCHP v11-02c). Geosci. Model Dev. 11, 2941–2953.

• Huynh, H. T. (1996). Schemes and constraints for advection. Fifteenth International Conference on Numerical Methods in Fluid Dynamics, Monterey, CA, Springer, 498–503.

• Krol, M., et al. (2005). The two-way nested global chemistry-transport zoom model TM5. Atmos. Chem. Phys. 5, 417–432.

• Lin, S.-J. (2004). A "vertically Lagrangian" finite-volume dynamical core for global models. Mon. Weather Rev. 132, 2293–2307.

• Lin, S.-J., and Rood, R. B. (1996). Multidimensional flux-form semi-Lagrangian transport schemes. Mon. Weather Rev. 124, 2046–2070.

• Putman, W. M., and Lin, S.-J. (2007). Finite-volume transport on various cubed-sphere grids. J. Comput. Phys. 227, 55–78.

• Russell, G. L., and Lerner, J. A. (1981). A new finite-differencing scheme for the tracer transport equation. J. Appl. Meteorol. 20, 1483–1498.