Tracer Conservation in Single Precision

AtmosTransport runs production transport in Float32. Data-centre GPUs such as the NVIDIA L40S execute Float64 at about 1/64 of the Float32 rate, Apple Metal has no Float64 kernels, and single precision halves the memory and bandwidth of the tracer, air-mass and flux fields. Float32 resolves about seven significant digits, and a year of transport applies tens of thousands of operator steps to every cell. This page describes how the model nevertheless conserves global tracer mass in Float32 to about 1e-6 of the burden per year, how that was measured, which numerical techniques it relies on, and what remains.

How single precision loses mass

Round-to-nearest is unbiased for a single operation. Global tracer mass drifts only when rounding errors share a sign across many cells and steps. Three mechanisms caused such drifts in the model:

1. Absorbed increments. Adding less than half an ulp of the target leaves the target unchanged, and adding a few ulps rounds the same way every time. A surface flux added to a 400 ppm CO₂ cell, or a weak diffusive exchange, loses or gains the same fraction on every step.

2. Near-identity linear solves. Implicit convection and diffusion solve (I + B) x = q with ‖B‖ ≪ 1. Elimination in Float32 rounds the diagonal 1 + bₖₖ and every update of a large qₖ by a small correction. Each solve then changes the column total by about 1e-8, in proportion to the tracer mass and with a fixed sign.

3. Geometry and constants rounded before use. A cell area, regridding weight, decay factor or time step computed in Float32 carries a fixed relative error. That error repeats on every application, for example in every area-weighted emission.

Flux-form advection applies each face flux with opposite signs to its two cells, including across cube-panel seams. The global total therefore changes only through the rounding of updated cell values. That rounding is mostly noise, with a small residual bias (see What remains).

How it was measured

Each run carries a uniform 400 ppm background tracer without sources. Its global mass is invariant, so any change is arithmetic error. SF₆ starts from its observed state. For tracers with emissions, a Float32 run is compared with a Float64 run of the same code. Global totals are the <tracer>_total_mass output variables, accumulated during the run in Float64 with compensated summation, independent of the storage precision. Runs switch operators on one at a time. For the noise analysis, the 3-hourly increments of a total are tested for a non-zero mean (t = mean divided by its standard error).

Archive	Levels, windows	Diffusion	Convection	Length
ERA5	C90 L66, hourly	TM5 dkg	TM5, collaborative LU	3 days, no sources
MERRA-2 (GEOS-Chem archive)	C90 L72, hourly	GEOS-Chem non-local VDIFF	CMFMC, DQRCU cloud base	7 days, with sources

Float32 runs used an NVIDIA L40S and Float64 runs an A100, starting 2021-12-01. In the Float64 MERRA-2 runs the background total stays constant to the last bit with every operator set. A Float32 run starts 2.5e-8 below its Float64 twin because Float32(400e-6) is 2.5e-8 low. That offset is constant and is removed from the comparisons below.

Results

ERA5, Float32, 3 days, no sources: relative change of global mass.

Operators	Background, before → after	SF₆, before → after
advection	+8.6e-9	−1.2e-9
advection + dkg diffusion	+9.0e-9 → +8.4e-9	−1.24e-8 → +2.5e-9
advection + TM5 convection	−9.7e-7 → +7.0e-9	−9.9e-7 → −2.0e-9
all operators	−9.7e-7 → +9.2e-9	−1.0e-6 → +2.1e-9

MERRA-2, 7 days with sources: Float32 − Float64. Background, CO₂ and SF₆ are relative to the global burden. Fossil CO₂ starts at zero and is relative to the emitted mass.

Operators	Background	CO₂	SF₆	Fossil CO₂
advection, before	+8.0e-9	−5.5e-9	−4.6e-9	+8.2e-7
advection, after	+8.0e-9	−4.2e-9	−4.1e-9	−4.0e-8
all operators, before	−9.6e-9	−2.4e-7	+1.1e-7	+8.5e-7
all operators, after	−2.3e-9	−2.0e-8	−3.7e-9	−1.9e-8

Over the week, the CO₂ burden changes by 3.0e-4 through net fluxes and the SF₆ burden by 6.9e-4. Before the fixes, Float32 therefore lost 8e-4 of the net CO₂ flux, almost all of it in the non-local boundary-layer scheme (advection + diffusion: −2.6e-7; advection alone: −5.5e-9).

Individual mechanisms, measured directly.

Quantity	Float32 before	After
C90 cell area, max. relative error	1.3e-4	5.7e-8 (½ ulp)
C90 Δx, Δy, max. relative error	1.8e-5	4.4e-8
Rn-222 decayed fraction per step (Δt = 300–900 s)	up to 3.4e-5	≤ 3.2e-8
GridFED → C90 regridding, F32 vs F64 global flux	+8.3e-7	6e-10 (rates rounded)
GridFED → C90 regridding, Float64 loss	2.1e-6	0
Model clock, 7 steps per hour, after one year	3.1 h behind	exact at window ends

A Float32 clock happens to stay exact over two years with 300, 400 or 900 s steps, but adaptive substepping can choose any step count; with seven steps per hour it is 24 h ahead after two years.

Techniques

Kahan-compensated accumulation

Every surface source keeps a persistent per-cell compensation c. The surface-layer update is

y = x − c;  t = s + y;  c = (t − s) − y;  s = t

so increments below the target's resolution carry forward instead of being lost (Kahan, 1965). After n steps the field has gained n·x + c, up to unbiased roundings of about eps·x per step; c is mass deposited ahead of schedule (negative: behind).

Error-free transformations

TwoSum (Møller, 1965; Knuth, 1969) returns s = fl(a + b) and the exact error e = (a + b) − s using six Float32 additions. GEOS-Chem's non-local boundary-layer scheme deposits fresh emissions over several layers. Each upper-layer addition is a TwoSum, so the mass that rounding kept out of a layer is known. The surface layer receives the emission minus what actually landed above, through the Kahan update. Each column thus gains the emitted mass even when the profile fractions do not sum to one in Float32.

On NVIDIA GPUs, ptxas fuses a multiply and a dependent add into one FMA (verified in the generated SASS for the deposit). The share x·fₖ may then enter the sums unrounded, and the bookkeeping is exact to about ½ ulp of the share per layer, with random sign.

Column mass ledger

The TM5 convection matrix and the dkg diffusion exchange conserve column mass exactly: their columns sum to one. Their Float32 solves do not. CMFMC convection's explicit update conserves too, since its interface fluxes cancel in the column sum, but it writes each layer's new mass q·m in Float32 (the cloud-base layer twice), and a change below half an ulp of the stored mass vanishes. Around each solve or update, the ledger

1. sums the participating cells before and after with compensated summation (Neumaier, 1974, in the branch-free TwoSum form of Ogita et al., 2005);

2. forms the residual r = (s₀ − s₁) + (c₀ − c₁). The solve changes the column total by far less than a factor of two, so s₀ − s₁ is exact (Sterbenz's lemma);

3. adds r and the column's carry c to the largest participating cell q with two error-free sums, and keeps what that cell cannot hold as the new carry:

t + t_lo = r + c,   q′ + s_lo = q + t,   c′ = fl(s_lo + t_lo).

The residual is a few ulps of the largest cell for TM5 convection (1–2e-8 of the active column) and about one ulp for dkg diffusion (≈5e-10 of the column). Spread over the column it would round away again. Added to the largest cell, it is rounded to that cell's spacing: the rounding error s_lo, at most half an ulp, is lost, and all of r is lost where it is below half an ulp, in 2–10 % of TM5 columns and 20–30 % of dkg columns. For a tracer peaked at the surface s_lo has one sign, so without the carry the losses added up: in one day of dkg diffusion alone, ERA5 C90 L117 fossil CO₂ lost 8.1e4 kg, 2.5e-9 of its burden (v0.5.0 and earlier).

With the carry, each cell is still rounded once per solve, and q′ + c′ = q + r + c up to one rounding of c′; with a zero carry the cell is bit-identical to the ledger without one. The carry is at most half an ulp of the cell plus the rounding of r + c; it is kept per column and tracer in the operator's workspace (TM5Workspace, CMFMCWorkspace, DiffusionWorkspace) and added with the column's next residual. The field plus the carries therefore conserves each column, and the field alone differs from it by the outstanding carries, which do not accumulate. Moving r to a smaller cell, where it is represented, also stops the drift but changes that cell by a larger fraction of its value, trading per-cell accuracy for column accuracy. In 200 repeated Float32 solves of surface-peaked test columns, the column totals drift by 7.7e-9 (TM5) and 1.9e-9 (dkg) without the carry and by 5e-11 and 8e-11 (the outstanding carries) with it; the largest error of a cell against Float64 stays 7.5e-6 and 3.4e-6 (test/core/test_float32_conservation.jl).

In one day of an ERA5 C90 L117 CATRINE run (2023-05-15, Float32 on an L40S, restarted from the v0.5 chain) with dkg diffusion only and no sources, the fossil-CO₂ total now changes by +57 kg (+1.7e-12) instead of −8.1e4 kg, and its 3-hourly values fluctuate within ±130 kg without a trend: the outstanding carries, about half an ulp of each column's largest cell at most. With TM5 convection only, the change is −300 kg instead of −780 kg, again without a trend. With all operators and sources, the fossil-CO₂ residual of the day falls from −2.7e-9 to −1.1e-10 of the burden, below that of advection and emissions alone (−2.7e-10). The carry costs one load, one store and a second TwoSum per column and tracer; on an L40S the time of a C90 L117 call with ten tracers is unchanged within noise (dkg 1.56–1.57 ms, TM5 convection 53–54 ms). The ledgers rely on IEEE arithmetic without fast-math reassociation, which would break TwoSum; CPU and CUDA are tested, while the Metal checks (test/diagnostic/test_metal_kernels.jl) still have to be run on a Mac.

Carries start at zero and are not written to restart files: a restarted run drops about half an ulp of each column's largest cell, once. A column that is not ledgered (no convection, no exchanging layer, or a zero-mass layer in a dkg column, which keeps its absorbing convention) keeps its carry for its next ledgered solve. Single-tracer dkg calls name their carry slot (tracer = idx, as the per-tracer Lin–Rood path does); without it they keep no carry.

The ledger corrects rounding, not physics. A residual larger than 16·n·eps·|q_max|, for n participating cells with largest value q_max, is left in place and never enters the carry. That bound is about 1e-6 of the column in Float32, far above the measured residuals, so a matrix that does not conserve mass still shows up in the budgets; a leak smaller than the bound would be returned like rounding, which is why the ledger runs only where the operator conserves by construction. A carry above the same bound is dropped; since a carry is at most about half an ulp of the cell it came from, that happens only when the column's largest cell shrank more than about 32n-fold since, for example when the column was emptied. A column overwritten with values of similar size keeps its carry, deposited once: a run restarted from reloaded fields into an existing model (reused workspaces) is therefore not bit-identical to one built fresh, by about half an ulp per column.

The ledger runs only where the operator conserves mass by construction. The linear-algebra routines stay plain solves. TM5 convection uses the cloud levels; with level merging (n_merge > 1) its ledger covers the fine layers that aggregation, coarse solve and split change, because the split also rounds every fine layer. dkg diffusion and the flux form of Kz diffusion use layers that exchange mass and leave isolated layers bit-exact. CMFMC convection uses the whole column where nothing leaves it: the cmfmc flux through the model top is exactly zero and every layer holds air (m and m/area above the operator's tiny); the clamped variant (clamp = true), which clips negative values and rescales the column, is not ledgered. The adjoint's forward replay applies the same ledgers without a carry: every replayed solve starts from a zero carry, so it is a function of the recorded state alone and checkpointed replays reproduce the recorded states bit for bit. The adjoint itself transposes the exact operator.

A delta-form solve, (I + B) Δ = −B q followed by q += Δ, would remove the bias at its source and also improve local accuracy. It needs the matrix before factorisation for every tracer batch, which exceeds the shared-memory budget of the collaborative GPU kernel.

Flux-form update of the Kz diffusion

Implicit diffusion with a Kz field solves the backward-Euler system A q′ = q for the new mixing ratio, with b_k = 1 + Δt (D_{k−½} + D_{k+½})/m_k on the diagonal and −Δt D_{k∓½}/m_k beside it (D the interface exchange in kg s⁻¹). Row k times m_k states m_k q′_k = m_k q_k + F_{k−½} − F_{k+½} with F_{k+½} = Δt D_{k+½} (q′_k − q′_{k+1}), so writing back m_k q′_k conserves the column in exact arithmetic. In Float32 the computed q′ solves the system with rounded entries. In the pure-pressure levels near the model top, Δt D/m_k is the same in every column (mass and exchange both scale with the cell area; exactly so with isothermal layer thickness and a horizontally uniform Kz, as in the golden configurations, approximately otherwise), so the rounding of b_k has the same sign in every column. The solve works on the departure from the column minimum, which is non-negative, so the error of the column total keeps one sign as well: it acts like a matrix whose columns do not sum to one. In the Float32 lat-lon L72 golden run it cost a tracer without background 2.0e-6 of its mass per day, through about 376 solves a day (finding OF-113).

The tracer mass is therefore updated from the fluxes, not from m_k q′_k: each F_{k+½} is computed once from the solved departures and applied with opposite signs to its two layers, rm′_k = (rm_k + F_{k−½}) − F_{k+½}. The column total telescopes; only the two roundings per layer change it. Those carry no common sign in general, but where the fluxes keep one direction, a flux between half an ulp of the receiving layer and half an ulp of the giving one is recorded by one layer only. Adding the net flux F_{k−½} − F_{k+½} in one step makes this worse: a small net change rounds into a large background the same way on every solve (CO₂ in the run below: +1.7e-7 per day instead of −9.4e-9). Even the two additions leave a small drift of one sign where the fluxes keep one direction (CO₂ in the table below, −9.4e-9 per day). The column ledger with its carry, as for dkg diffusion, returns it. The solve and the fluxes are those of the matrix form, the two agree in exact arithmetic, and the adjoint transposes the same operator. Layers without exchange keep their mass bit for bit.

The flux form trades local accuracy for column conservation. In floating point its result differs from the matrix form's by m_k r_k, with r the residual of the computed tridiagonal solve, which grows like ε (1 + 2x_k) with x_k = Δt (D_{k−½} + D_{k+½})/m_k: the matrix form puts its rounding into the column total, the flux form into the individual layers. In a Float32 probe of strongly mixed boundary-layer columns (x up to about 1e3) against an exact backward-Euler solve, the worst relative error of a single cell was 1.4e-4 for the flux form and 1.0e-5 for the matrix form, while the column mass error per solve was 3e-8 and 5e-6. The next strongly mixing solve damps this per-cell noise, so it does not accumulate, and cells stay non-negative unless x approaches 1/ε. It is also why the golden fields of the affected Float32 cases change by 2.5e-5 to 8.3e-5 of their maximum, more than a few ulps. In Float64 the two forms differ by about 5e-16 for x ≲ 1 and 5e-14 for x ≈ 1e3.

Diffusion's own change of each global total in one day of the lat-lon L72 golden run (run_ll72_upwind_cpu_f32, 376 solves), relative to the burden (fossil CO₂: to the emitted mass):

Tracer	matrix form	flux form	flux form, ledger without carry	flux form, ledger and carry
blob (no background)	−2.0e-6	+1.2e-8	+2.7e-10	−1.5e-10
CO₂ (with surface fluxes)	−2.6e-9	−9.4e-9	+2.7e-10	+1.5e-11
fossil CO₂	−2.6e-7	+4.2e-9	+2.4e-8	−3.5e-10
Rn-222	−7.9e-8	+1.3e-8	+1.7e-9	+4.6e-10
uniform 400 ppm	+2.0e-9	−1.0e-10	−1.0e-10	−4.3e-11

The flux form alone removes the bias of the solve but leaves about 1e-8 per day, more than the matrix form for CO₂; without the carry the ledger leaves a fossil-CO₂ drift, as it did for dkg diffusion. With both, diffusion changes each total by less than 5e-10 per day. In repeated Float32 solves of test columns with pure-pressure upper levels, the field plus the outstanding carries conserves the global total to 2e-15, against 7.7e-7 to 5.1e-6 for the matrix form in 192 solves (test/core/test_diffusion_flux_form.jl, CPU and CUDA).

Float64 geometry and constants, rounded once

Quantities computed once and used many times are evaluated in Float64 on the host and rounded to the run precision a single time:

• cubed-sphere corners, cell areas and edge lengths;

• conservative-regridding geometry. Meshes of either precision regrid on a Float64 sphere and share cached weights. Source grids stored in Float32 (GridFED) take their spacing from the full span, and extents within 1e-3 of a cell of ±90° or 360° snap to them. Otherwise the regridder misses a sliver of the sphere;

• the decay decrement d = expm1(−λΔt), applied in kernels as c + c·d;

• the model clock, computed from the window and step counters, so window ends fall exactly on multiples of the window length.

Kernels never see Float64 values. They stay type-stable in the run precision, as Metal requires and the L40S rewards.

Preprocessing builds the same mesh. Binaries written with float_type = "Float32" before this change used Float32 cell areas, with errors up to 1.3e-4 per cell at C90, 6.8e-4 at C180 and 9.7e-3 at C720 (global area −1.6e-8, −3.1e-7 and +2.1e-6), wherever the preprocessor converts between pressure and mass. Such binaries should be regenerated, especially at C180 and finer.

Diagnostics in Float64

Totals meant to reveal 1e-8 drifts cannot be Float32 reductions. Snapshot totals, total_mass and the log totals share one compensated Float64 reduction: in place on CPU and CUDA, in bounded host slabs on Metal.

What remains

Advection. The 3-hourly change of the background total has a standard deviation of 5.5e-10. Its mean is +1.4e-10 (MERRA-2, t = 1.9) to +3.6e-10 (ERA5, t = 3.2), i.e. at most about +1e-6 per year. SF₆ shows no detectable mean (t = −1.0). The cause of the bias for the uniform background has not been identified. Removing storage rounding altogether would need double-single tracer storage (two Float32 words per value), doubling memory traffic. That is not warranted at this level.

CMFMC convection. In the MERRA-2 week above, the remaining CO₂ difference with all operators (−2.0e-8 of the burden, about −1e-6 per year) matched what CMFMC convection added to advection alone (−2.2e-8); CMFMC had no column ledger then (v0.5.0 and earlier). In one day of the MERRA-2 C90 L72 CATRINE run (2023-05-15, Float32 on an L40S, restarted from the v0.5 chain, no sources), CMFMC convection alone changed the SF₆ total by +1.3e-9 of its burden and background CO₂ by +2.5e-9 without the ledger, and by −5.5e-12 and +4.9e-12 with it. With both ledgers, every operator and the sources, the day's SF₆ residual falls from +3.8e-9 to −6.6e-10 of the burden, the level of advection alone (−6.1e-10). In 200 repeated updates of a surface-peaked test column, the ledger also keeps the cells closer to Float64: within 1.3e-6 instead of 1.7e-5, and a uniform mixing ratio within 4.3 instead of 76 eps (relative; test/core/test_float32_conservation.jl). A C90 L72 CMFMC call with ten tracers takes about 1 % longer on an L40S (5.42 → 5.47 ms).

Emission rates. Rounding each cell's rate to Float32 changes global emission totals by less than 1e-9.

Initial conditions. Converting mixing ratios to tracer mass rounds once (the 400 ppm offset above). Output fields are stored in Float32; budgets must use the _total_mass variables, not sums of output fields.

Not covered. The ledger does not cover the clamped CMFMC variant, at about 5e-6 per month, which the production configurations do not use. Kz diffusion (local-Kz fields, and the lat-lon and reduced-Gaussian runs) is covered since its flux-form update (see above).

Global mass fixers, which rescale tracer fields to a target total, would hide operator errors instead of removing them and are not used.

Reproducing

The experiment is scripted in scripts/diagnostics/float32_conservation_attribution.jl:

julia --project=. scripts/diagnostics/float32_conservation_attribution.jl write-configs --out DIR
julia --project=. scripts/run_transport.jl DIR/merra2_full_float32.toml    # each of the 16 configs, on a GPU
julia --project=. scripts/diagnostics/float32_conservation_attribution.jl budget --out DIR

write-configs builds the 16 run configurations (two archives, four operator sets, two precisions), checks each with validate_config and writes DIR/<tag>.toml; the inputs are taken from --data-root (default ~/data/AtmosTransport). budget reads the <tracer>_total_mass series of every run found and prints the tables above: the Float32 drift, Float32 − Float64 at the last common snapshot net of the initial offset (fossil CO₂ relative to the emitted mass), and the mean, standard deviation and t of the 3-hourly increments of the source-free tracers (twin_relative_difference, drift_statistics). --tag and --prefix name the run folders and snapshot files: a template with the placeholders {archive}, {set}, {float_type} and {ft}, for all runs, for one archive (merra2=...) or for one archive and precision (merra2:float32=...). The runs behind this page (memo docs/memos/2026-10-06_float32_conservation_attribution.md; the "after" runs are its era5_*_final and full_final2, the latter compared with the Float64 run full_final) are read with

julia --project=. scripts/diagnostics/float32_conservation_attribution.jl budget \
    --out ~/data/AtmosTransport/output/precision_attribution --prefix precision \
    --tag 'merra2:float32={set}_final2_f32' --tag 'merra2:float64={set}_final_f64' \
    --tag 'era5=era5_{set}_final_{ft}'

which prints the "after" values of the ERA5 table and the row "all operators, after" of the MERRA-2 table. merra2={set}_regridfix_{ft} gives the row "advection, after", and the templates {set}_{ft} and era5_{set}_{ft} give the values "before".

Longer twin runs are compared with scripts/diagnostics/compare_float32_float64_runs.jl: per tracer, the relative difference of the global totals at every snapshot and, once a day, the area-weighted mean and RMS of the column-mean difference (weighted_difference_stats). It writes float32_vs_float64.json and float32_vs_float64.nc, from which scripts/visualization/plot_float32_vs_float64.py draws the figures. For the MERRA-2 C90 CATRINE configuration over 761 days (December 2021 to December 2023), Float32 − Float64 of the global totals, net of the initial offset, is +9.0e-7 for CO₂, +1.6e-7 for fossil CO₂ (relative to the emitted mass), +1.1e-6 for SF₆ and −3.6e-8 for Rn-222; the column-mean RMS difference on the last day is 0.015 ppm for CO₂ and 0.0009 ppt for SF₆.

test/core/test_float32_conservation.jl checks the error-free and compensated sums, that the ledger leaves a leaking matrix visible, the ledger carry (its deposit, repeated-solve drift of TM5 and dkg columns against Float64, and that the TM5 workspaces keep it on every topology), the geometry, the decay decrement, the source-grid extents and the profile deposit. test/core/test_tm5_hessenberg.jl and test/helpers/conservative_dkg.jl require the convection and diffusion ledgers to conserve each column to within one ulp of its largest cell. test/core/test_multiday_source_clock.jl checks that the clock lands on window ends with seven steps per window.

References

• Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, ch. 4.

• Kahan, W. (1965). Further remarks on reducing truncation errors. Commun. ACM 8, 40.

• Knuth, D. E. (1969). The Art of Computer Programming, Vol. 2, Seminumerical Algorithms, §4.2.2.

• Møller, O. (1965). Quasi double-precision in floating point addition. BIT 5, 37–50.

• Neumaier, A. (1974). Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen. Z. Angew. Math. Mech. 54, 39–51.

• Ogita, T., Rump, S. M., and Oishi, S. (2005). Accurate sum and dot product. SIAM J. Sci. Comput. 26, 1955–1988.