Analysis API

Analysis holds the shared numerics of post-run diagnostics: reading a time series from daily snapshot files, area-weighted means on equal-area latitude bands, the seasonal, trend and interannual decompositions used to compare transported tracers with observations, skill statistics, and column profile tools. Diagnostic command lines under scripts/diagnostics/ take these numerics from here, and plotting scripts read the files those command lines write. The names are exported from AtmosTransport.Analysis (using AtmosTransport.Analysis), not from the top level.

Conventions shared by every function:

• arrays are read in NCDatasets order, for example (Xdim, Ydim, nf, time) for a cubed-sphere column mean, and time is the last dimension of a series;

• sums run in Float64 in a fixed (linear-index) order, so a result does not depend on the input's element type or on threading;

• NaN marks a missing value: band means and skill statistics leave it out, while the decompositions propagate it (and seasonal_trend_decomposition rejects it at observed months);

• statistics take the model first and the reference second and use population (1/n) moments;

• column profiles are ordered top first (layer 1 at the model top), as in AtmosTransport snapshots.

A typical Hovmöller comparison reads the daily column means with daily_snapshot_series, forms (nbands, ntime) band means with one BandAssignment, averages them to months, and passes the stacked model and observation series to seasonal_trend_decomposition and skill_metrics.

The global totals a run writes at every snapshot (<tracer>_total_mass, a compensated Float64 sum) feed two conservation diagnostics. drift_statistics measures the change of one run's total against its initial value and tests its per-snapshot increments for a mean ($t$); for a tracer without sources this is pure arithmetic error. twin_relative_difference compares a run with its twin in another precision frame by frame, net of the rounding of the initial state. weighted_difference_stats gives the area-weighted mean and RMS of a field difference, such as the column-mean difference of the two runs. scripts/diagnostics/float32_conservation_attribution.jl and scripts/diagnostics/compare_float32_float64_runs.jl are built on them; see Tracer Conservation in Single Precision.

GEOSChemCATRINEArchive reads the output of the GEOS-Chem (GCHP) CATRINE reference run, one GEOSChem.CATRINE_inst.YYYYMMDD_HHMMz.nc4 file per 3-hourly output time, flat or under YYYY/MM/. global_totals turns a file into global sums in GEOS-Chem's dry-air-equivalent storage units: the dry-air mass $\sum \mathrm{Met}\mathrm{\_}\mathrm{AD}$, the burden $\sum \chi \mathrm{Met}\mathrm{\_}\mathrm{AD}$ of each tracer and its emission rate $(\sum E\mathrm{Met}\mathrm{\_}\mathrm{AREAM}2)M_{\mathrm{air}}/M$, with GEOS-Chem's molar masses (Parameters.GEOSCHEM_SPECIES_MOLAR_MASS, Parameters.DRY_AIR_MOLAR_MASS). A missing, truncated or incomplete file reads as nothing and is listed by unreadable_files. Files are read into memory on all threads and decoded under the package's NetCDF lock. scripts/diagnostics/catrine_true_mass_balance.jl --model geoschem builds GEOS-Chem's true mass balance on it.

budget_closure checks the global budget of every tracer of a finished run. The runtime conserves tracer storage $S=\sum _c{\chi }_c{m}_c$ (dry mixing ratio times carrier air mass, the unit of <tracer>_total_mass), and transport drops out of its global sum, so between the first and last snapshot

$$\mathrm{\Delta }S=S(t_1)-S(t_0)=E(t_0,t_1),$$

with $E$ the storage the tracer's surface sources added. $E$ comes from the run's own sources (RuntimeSurfaceSources: built by the runner's loader on the run's grid and integrated with Operators.integrated_emission), so cell areas, unit factors, the storage scale $M_{\mathrm{air}}/M$ and slice times are the run's and molar masses cancel. A tracer passes when $|\mathrm{\Delta }S-E|\le \rho |E|+k{\epsilon }_{32}({\chi }_0)M$, where ${\epsilon }_{32}({\chi }_0)$ is the Float32 spacing of its initial carrier and $M$ the run's initial air mass (BudgetTolerance: strict $\rho ={10}^{-4}$, $k=1$; operational $\rho ={10}^{-3}$, $k=4$), its values are finite and, with a positive carrier, no column mean drops to zero. RunContext holds the run as the runner builds it from its config. carrier_drawdown gives, before a run, the largest cumulative loss of global mixing ratio each tracer's sources cause, to compare with its carrier. scripts/diagnostics/tracer_budget_closure.jl series <run.toml> writes the JSON report (report_dict) and exits 1 when the gated level fails. Another way to compute $E$ plugs in as a subtype of AbstractExpectedEmission.

AtmosTransport.Analysis Module

Analysis (src)

Shared numerics for post-run diagnostics: reading daily snapshot series, equal-area latitude-band means, seasonal and interannual decompositions, skill statistics, column profile tools (column-mass band weights, log-pressure interpolation), the drift and twin-run differences of global tracer totals, a reader for the output of the GEOS-Chem CATRINE reference run (GEOSChemCATRINEArchive), and the global budget closure of a finished run (budget_closure, carrier_drawdown).

Diagnostic command lines in scripts/diagnostics/ take these numerics from here instead of carrying their own copies, and plotting scripts read the NetCDF or JSON those diagnostics write. Every function works on plain arrays on the host and accumulates in Float64 in a fixed order, so a result does not depend on threading or on the element type of the input.

src/AtmosTransport.jl loads this module after Models and before Downloads, so that diagnostics can build on the runtime's own types; no runtime module depends on it. The functions here use NCDatasets, Dates, LinearAlgebra and Statistics, the compensated sum of Architectures, the constants of Parameters and the NetCDF lock of Output; the budget closure builds its view of a run (RunContext) and the expected emissions with the runner's own config parsing, grids, output naming and surface-flux sources. See src/Analysis/README.md for the file map and docs/src/api/analysis.md for the reference.

source

AtmosTransport.Analysis.CATRINE_GEOSCHEM_TRACERS Constant

CATRINE_GEOSCHEM_TRACERS

The four tracers of the GEOS-Chem CATRINE reference run, labelled as in the AtmosTransport CATRINE runs: co2_natural (SpeciesConcVV_CO2, EmisCO2_Total), co2_fossil (SpeciesConcVV_FossilCO2, Emis_FossilCO2_Total), sf6 (SpeciesConcVV_SF6, EmisSF6) and rn222 (SpeciesConcVV_Rn222, EmisRn_Soil, half-life RN222_HALF_LIFE). Molar masses are GEOS-Chem's (GEOSCHEM_SPECIES_MOLAR_MASS; the CO₂ value for both CO₂ tracers).

source

AtmosTransport.Analysis.AbstractExpectedEmission Type

AbstractExpectedEmission

How budget_closure obtains the expected storage change of a tracer. A subtype implements expected_storage_change(method, ctx, spec, t0, t1). The default, RuntimeSurfaceSources, integrates the run's own surface-flux sources; a validation harness can plug in another method, for example a reimplementation of a retired script with its constants.

source

AtmosTransport.Analysis.AbstractLatitudeBands Type

AbstractLatitudeBands

A partition of latitudes into bands. A subtype implements nbands, band_edges and band_index(bands, lat), which returns the band of a cell centre latitude in 1:nbands(bands), or 0 for a latitude in no band (such a cell enters no mean).

source

AtmosTransport.Analysis.BandAssignment Type

BandAssignment(bands, lat, area)

The band of every cell of a grid and the cell areas, built once and reused by band_means. lat (cell centre latitudes, degrees north) and area (cell areas, or any non-negative weights) have the shape of the grid, for example (Xdim, Ydim, nf) for a cubed sphere read with NCDatasets or (Nx, Ny) for a lat-lon grid. Areas must be finite and non-negative.

source

AtmosTransport.Analysis.BudgetReport Type

BudgetReport

Result of budget_closure:

• config_path (empty when the config was passed as a Dict) and files, the snapshot files read, in time order;

• start and stop: the times of the first and last snapshot, between which the budgets are formed; frames: the number of distinct snapshot times;

• tolerance (BudgetTolerance) and emission, the method of the expected changes (repr of an AbstractExpectedEmission);

• tracers: one TracerBudget per checked tracer, in the order of the run's tracer specifications;

• skipped: tracers not checked, with the reason (written to no snapshot, or a chemistry sink that snapshots do not resolve);

• passed and operational_passed: every checked tracer passes that level; max_relative_error: the largest relative_error.

report_dict turns it into the JSON report of scripts/diagnostics/tracer_budget_closure.jl.

source

AtmosTransport.Analysis.BudgetTolerance Type

BudgetTolerance(; strict = 1e-4, operational = 1e-3,
                  strict_carrier_ulps = 1, operational_carrier_ulps = 4)

Tolerances of budget_closure. A tracer with expected storage change $E$, actual change $\mathrm{\Delta }S$, initial carrier ${\chi }_0$ and air mass $M$ passes a level with relative tolerance $\rho$ and $k$ carrier ulps when

$$|\mathrm{\Delta }S-E|\le \rho |E|+k{\epsilon }_{FT}({\chi }_0)M,{\epsilon }_{FT}(x)=\mathrm{eps}(\mathrm{FT}(x)).$$

The strict level (${10}^{-4}$ plus one carrier ulp) is the check of the retired check_trendy_npp_rh_run.py, which gated the NPP/Rh campaign (scripts/run_trendy_npp_rh_campaign.py still gates on it); the operational level (${10}^{-3}$ plus four ulps) is the guard of scripts/run_transcom_campaign.py. Every value must be finite and non-negative.

source

AtmosTransport.Analysis.EqualAreaBands Type

EqualAreaBands(n)

n latitude bands of equal area. On a sphere the area south of latitude φ is proportional to $1+\sin \phi$, so equal areas mean edges equally spaced in $\sin \phi$:

$$s_k=\frac{2k}{n}-1,k=0,\dots ,n,$$

evaluated as (k / n) * 2 - 1 in Float64 (the rounding of the Python diagnostics this replaces). Band $b\in 1:n$ holds the latitudes with $s_{b-1}\le \sin \phi <s_b$, where $\sin \phi$ is evaluated as sin(deg2rad(φ)) and compared in IEEE arithmetic (so -0.0 lies with 0.0, as in NumPy digitize); the north pole ($\sin \phi =1$) belongs to band n. Every latitude in $[-90°,90°]$ is in exactly one band.

source

AtmosTransport.Analysis.GEOSChemCATRINEArchive Type

GEOSChemCATRINEArchive(dir; tracers = CATRINE_GEOSCHEM_TRACERS, interval = Hour(3))

The output of the GEOS-Chem CATRINE reference run under dir: one file GEOSChem.CATRINE_inst.YYYYMMDD_HHMMz.nc4 per output time, every interval, flat in dir or in dir/YYYY/MM/. tracers is a tuple of GEOSChemTracers.

The archive records the files it could not read: with_gc_file (and global_totals, which uses it) returns nothing for a missing, truncated or otherwise unreadable file or one that lacks a variable, and unreadable_files lists them with the reason. Reads are safe from several threads.

source

AtmosTransport.Analysis.GEOSChemTracer Type

GEOSChemTracer(name, concentration, emission, molar_mass; half_life = Inf)

One tracer of a GEOS-Chem run and the output variables that describe it:

• name: the diagnostic's label for the tracer (a Symbol, for example :co2_natural);

• concentration: the dry mole fraction variable [mol mol⁻¹], for example "SpeciesConcVV_CO2";

• emission: the surface emission flux variable [kg species m⁻² s⁻¹], for example "EmisCO2_Total";

• molar_mass: GEOS-Chem's molar mass of the species [kg mol⁻¹], from Parameters.GEOSCHEM_SPECIES_MOLAR_MASS;

• half_life: radioactive half-life [s]; Inf for a stable tracer.

source

AtmosTransport.Analysis.RunContext Type

RunContext(config_path::AbstractString)
RunContext(config::AbstractDict; config_path = "")

The parts of a transport run that post-run diagnostics need, built with the runtime's own code from the run config (a TOML path or the parsed Dict):

• config: the parsed run config;

• grid: the grid of the first transport binary, opened on the CPU in the run's float type ([numerics].float_type), as the runner builds it;

• binary_paths: the transport binaries, in run order ([input]);

• tracer_specs: the tracer specifications the runner parses ([tracers.<name>]), sorted by name; without [tracers], the legacy single tracer [run].tracer_name (default CO2) of a lat-lon or reduced-Gaussian run, initialised from [init] and without a surface flux;

• reference_time: the run start, [input].start_date at 00 UTC. Snapshot time values are hours since it, and flux slice times are seconds since it;

• windows: every met window the run steps through, in order, as (; binary, start, length, steps, step): the index of its binary, its start in seconds since the run start, its length (the header's window length in the run's float type), its number of advection steps (the binary's steps_per_window schedule) and the step length the runtime applies, FT(length) / FT(steps) in the run's float type (in a Float32 run seven steps of a 3600 s window add up to 3599.99994 s). The windows are those of [run].start_window through [run].stop_window of every binary, as the runner selects them;

• run_seconds: the run length, computed as the runner computes it from every binary's window layout;

• air_mass: the total air mass of window 1 of the first binary, in kg (dry air for the default dry-basis binaries): the air mass the runner builds the initial tracer state on, also when [run].start_window is larger;

• decaying_tracers: tracers with a [chemistry] sink (half-lives of kind = "decay").

floattype(ctx) returns the run's float type. Building a context reads the header of every binary (about 40 ms each) and the first window of the first. The run must have [input].start_date, and its window range and tracers must be ones the runner accepts (1 ≤ start_window ≤ stop_window ≤ the window count of each binary; start_window = 1 and every window of every binary for a run of several binaries; on a cubed-sphere grid start_window = 1 and a [tracers] table).

source

AtmosTransport.Analysis.RuntimeSurfaceSources Type

RuntimeSurfaceSources()

Expected storage change from the run's own surface-flux sources: the tracer's source is built by build_surface_flux_sources on the run's grid, in the run's float type, with the run start as reference time (so with the run's cell areas, unit factors, storage scale $M_{air}/M_{tracer}$ and slice times), and integrated over the run's advection steps (ctx.windows) with Operators.integrated_emission and Operators.emission_segments, so $[t_0,t_1]$ must consist of whole met windows, as the snapshot times do. The runtime starts step $j$ of a window of length $L$ and $n$ steps at $t_w+Lj/n$ and applies the source for a step $\mathrm{\Delta }t$ = FT(L) / FT(n) in the run's float type. Static, stepwise and conservative-mean sources are integrated exactly piece by piece, so over windows whose steps tile them ($n\mathrm{\Delta }t=L$) one exact integral gives the step sum; in the other windows (a Float32 run with seven steps per window, for example) the steps are summed one by one. A linearly interpolated source is sampled at every step start: window by window with integrated_emission and the window's step where the steps tile it (its starts $t_w+j\mathrm{\Delta }t$ equal the runtime's to Float64 rounding), step by step elsewhere. The source holds the slices the runner loads for the interval (through 1 s past $t_1$). Piece durations are Float64; the runtime rounds the duration $\mathrm{\Delta }t\cdot f$ of a piece of a step split at a slice time to the float type (dt_FT * FT(frac), with its blend weights in the float type too), a difference of a few float-type roundings of that piece's emission, the order of the rounding of each per-cell increment the operator adds; it is not reproduced here.

source

AtmosTransport.Analysis.TracerBudget Type

TracerBudget

Budget of one tracer in a BudgetReport. Masses are tracer storage in kg (mixing ratio × carrier air mass, the unit of <tracer>_total_mass).

• tracer, carrier (initial carrier mixing ratio ${\chi }_0$: the background of a uniform initial condition, otherwise the initial mean $S(t_0)/M$), air_mass ($M$, RunContext);

• initial_storage $S(t_0)$ and final_storage $S(t_1)$, the first and last snapshot totals;

• expected_change $E$ and actual_change $\mathrm{\Delta }S=S(t_1)-S(t_0)$;

• relative_error $=|\mathrm{\Delta }S-E|/max(|E|,1\mathrm{kg})$ and error_vmr $=(\mathrm{\Delta }S-E)/M$, the error as a mixing ratio;

• carrier_ulp $={\epsilon }_{FT}({\chi }_0)M$ for the run's float type, and the absolute strict_tolerance and operational_tolerance of BudgetTolerance;

• minimum_column_mean: the smallest column-mean mixing ratio in any snapshot (NaN when the files have no <tracer>_column_mean), column_checked: whether the files have it, and finite: whether every total and column mean is finite;

• strict_passed, operational_passed: the closure holds at that level, the values are finite and, for a tracer with a positive carrier, every column mean stays positive (the carrier is not exhausted). A tracer with a positive carrier whose column means were not written cannot show that and fails.

source

AtmosTransport.Analysis.band_edges Method

band_edges(bands::AbstractLatitudeBands) -> Vector{Float64}

The nbands(bands) + 1 band edges, ascending, in the coordinate the band set is defined in: $\sin \phi$ for EqualAreaBands.

source

AtmosTransport.Analysis.band_index Method

band_index(bands::AbstractLatitudeBands, lat) -> Int

Band of the latitude lat (degrees north) in 1:nbands(bands), or 0 when the latitude is in no band. A latitude outside $[-90°,90°]$ (or NaN) is an ArgumentError.

source

AtmosTransport.Analysis.band_means! Method

band_means!(out, field, assignment::BandAssignment) -> out

In-place band_means: out has size (nbands, trailing...), where trailing are the dimensions of field after the grid dimensions.

source

AtmosTransport.Analysis.band_means Method

band_means(field, lat, area, bands) -> Array{Float64}
band_means(field, assignment::BandAssignment) -> Array{Float64}

Area-weighted mean of field in each latitude band,

$${\overline{x}}_b=\frac{\sum _{c\in b,x_c\ne \mathrm{NaN}}{a}_c{x}_c}{\sum _{c\in b,x_c\ne \mathrm{NaN}}{a}_c},$$

with cell areas $a_c$. NaN marks a missing value: such a cell is left out of the numerator and the denominator, so a band whose cells are all NaN has the mean NaN. Infinite values are not missing and propagate.

The leading dimensions of field must have the shape of the grid (the shape of lat); trailing dimensions (time, series) are looped, and the result has size (nbands, trailing...): a Vector for a single field. Sums run in Float64 over the cells in linear (column-major) order, the order of the Python numpy.bincount diagnostics this replaces, so the same inputs give the same bits. Build a BandAssignment once when averaging many fields on one grid.

source

AtmosTransport.Analysis.budget_closure Method

budget_closure(run; output_glob = nothing, tolerance = BudgetTolerance(),
               emission = RuntimeSurfaceSources(), tracers = nothing,
               require_complete = true) -> BudgetReport

Check the global budget of every tracer of a finished run: the change of its storage between the first and last snapshot against the storage its sources added over the same interval.

run is the run config (a TOML path or the parsed Dict, as passed to run_driven_simulation) or a RunContext. For each tracer:

• the actual change $\mathrm{\Delta }S=S(t_1)-S(t_0)$ comes from <tracer>_total_mass in the run's snapshot files, read with daily_snapshot_series from the files the run names in [output].path (run_snapshot_paths), or from the files matching output_glob (sorted), for a run whose output was moved;

• the expected change $E$ is expected_storage_change(emission, ctx, spec, t0, t1);

• the tolerance is $\rho |E|+k{\epsilon }_{FT}({\chi }_0)M$ (BudgetTolerance), with $FT$ the run's float type, ${\chi }_0$ the tracer's initial carrier and $M$ the run's initial air mass;

• every total and <tracer>_column_mean value must be finite, and a tracer with a positive carrier must keep every column mean positive (so its column means must be in the snapshots, as in check_trendy_npp_rh_run.py).

With require_complete = true the snapshot times must be exactly those the run writes (the [output] schedule through the end of the run, hour 0 included when scheduled), so a missing file or an unfinished run is an error; with false the budget covers whatever frames the files hold. tracers restricts the check to some tracer names. Tracers missing from the snapshots and tracers with a [chemistry] sink (decay) are listed in skipped: a daily total does not resolve the decay integral.

The budget is in tracer storage, so molar masses cancel; a residual is a transport or accumulation error, not a unit error. Replaces check_trendy_npp_rh_run.py, whose expected change used the flux file's cell_area, M_air = 0.02896546 and whole-day sums; this one uses the runtime's sources, so it is valid for any run and any temporal scheme.

source

AtmosTransport.Analysis.carrier_drawdown Method

carrier_drawdown(run) -> Dict{Symbol, Float64}

Largest cumulative loss of global mixing ratio that each tracer's sources cause over the run, for checking before a run that a carrier (the initial background) is large enough to absorb a net sink such as NPP uptake:

$$D=max(0,-\underset{t_0\le t\le T}{min}\frac{E(t_0,t)}{M}),$$

with $E(t_0,t)$ the storage the tracer's runtime source adds from the run start $t_0=0$ to $t$ (RuntimeSurfaceSources), $T$ the run length and $M$ the run's initial air mass (the air mass of the first window of the first binary, RunContext). $E$ is accumulated piece by piece of Operators.emission_segments: at every slice time for static and stepwise sources, whose rate is constant between slice times, and in addition at every advection step for conservative-mean and linearly interpolated sources and in windows whose float-type steps do not tile them (the steps of RuntimeSurfaceSources). The minimum is taken at the ends of the pieces, so it is the run's minimum at its step ends, or lower where a slice time falls inside a step. A carrier ${\chi }_0$ covers the run when $D<{\chi }_0$. Tracers without a surface flux have $D=0$. run is a run config (path or Dict) or a RunContext; the run needs no output.

Replaces the npp_max_global_drawdown_ppm audit of prepare_trendy_s3_npp_rh.py, which used whole-day sums of the flux file, M_air = 0.02896546 and a fixed dry-air mass of 5.1353e18 kg.

source

AtmosTransport.Analysis.column_band_weights! Method

column_band_weights!(w, m, bands) -> w

In-place column_band_weights: w is length(m) × length(bands).

source

AtmosTransport.Analysis.column_band_weights Method

column_band_weights(m, bands) -> Matrix{Float64}

Weights that average a column profile over bands of column air mass.

m holds the air mass of each layer of one column, top first, in any mass unit (kg m⁻² or kg per cell); it must be finite and non-negative with a positive total. bands is a collection of (lower, upper) fractions of the column air mass counted from the surface, $0\le$ lower $<$ upper $\le 1$: (0, 0.1) is the lowest 10 % of the air mass, (0.9, 1) the top 10 % and (0, 1) the whole column. A NamedTuple of pairs works and keeps its order. The result is a length(m) × length(bands) matrix.

With the column mass $M=\sum _k{m}_k$ and the fraction of the column below and at the top of layer $k$,

$$f_k^{\mathrm{b}\mathrm{o}\mathrm{t}}=\frac{1}{M}\sum _{j>k}{m}_j,f_k^{\mathrm{t}\mathrm{o}\mathrm{p}}=f_k^{\mathrm{b}\mathrm{o}\mathrm{t}}+\frac{m_k}{M},$$

the weight of layer $k$ in band $(a,b)$ is the air mass of the part of the layer inside the band,

$$w_k=max(0,min(b,f_k^{\mathrm{t}\mathrm{o}\mathrm{p}})-max(a,f_k^{\mathrm{b}\mathrm{o}\mathrm{t}}))M,$$

so that $\sum _k{w}_k=(b-a)M$ and the mass-weighted band mean of a mixing ratio is $\sum _k{w}_k{q}_k/\sum _k{w}_k$. $M$ is summed top to bottom and $\sum _{j>k}{m}_j$ is formed as the running sum from the surface minus $m_k$, the arithmetic of the Python benchmark this replaces.

source

AtmosTransport.Analysis.daily_snapshot_series Method

daily_snapshot_series(paths, variables; origin=nothing, keep_frame = nothing)
    -> (; times, values)

Concatenate the frames of runtime snapshot files into one series per variable.

• paths are read in the order given; pass them in time order (the daily files of a run are named <prefix>_YYYYMMDD.nc, so sort suffices).

• variables names snapshot variables whose last dimension is time, for example "co2_natural_column_mean" (Xdim, Ydim, nf, time) or "air_mass_per_area" (Xdim, Ydim, nf, lev, time).

• origin is optional for current files: their CF time coordinate records hours since <run start> and is decoded directly. Pass origin only to read legacy files whose time units are the bare string "hours"; when supplied for a current file it must agree with the encoded origin.

Each frame is stamped origin + hours. A frame whose hour was already read is a duplicated boundary frame and is skipped: the first occurrence wins. A frame earlier than the latest frame read that is not such a duplicate means the paths are out of order or come from different runs, and is an error. The first file of a run holds the initial state at hour 0 as well as its end-of-period frame; both are kept. Callers that want end-of-day states only drop the frame at origin.

keep_frame, a function of the frame time returning a Bool, reads only the frames it accepts, for example t -> Dates.hour(t) == 0 for one field per day from 3-hourly snapshots. Duplicates and the time order are checked on every frame, kept or not. Without it every frame is read.

Returns times::Vector{DateTime}, strictly increasing, and values::Dict{String, Array}: for each variable an array with the variable's non-time dimensions in NCDatasets order (for example (Xdim, Ydim, nf)) followed by time. The element type is the stored floating-point type; values equal to _FillValue or missing_value become NaN. Every file must hold each variable with the same non-time shape and element type.

source

AtmosTransport.Analysis.doy_linear_detrend Method

doy_linear_detrend(x, dates) -> Array{Float64}

Interannual variability of daily series: for each calendar day, fit and remove a straight line across years.

dates holds one date per entry of the last dimension of x (time); all other dimensions of x are independent series. Samples are grouped by calendar day (month and day of month), with 29 February folded onto 28 February, so in a leap year that key has two samples of the same year. For the samples $i\in I_d$ of day $d$, in years $y_i$,

$$x_i\approx {\alpha }_d+{\beta }_d(y_i-{\overline{y}}_d),r_i=x_i-{\alpha }_d-{\beta }_d(y_i-{\overline{y}}_d),$$

with ${\overline{y}}_d$ the mean year of the samples; the least-squares fit gives ${\alpha }_d$ the mean of the samples and ${\beta }_d=\sum (y_i-{\overline{y}}_d)x_i/\sum (y_i-{\overline{y}}_d{)}^2$. The result is the residual $r$. When all samples of a day fall in one year the minimum-norm solution has ${\beta }_d=0$, and a day seen once has residual 0.

This removes the mean seasonal cycle and any linear drift of it over the record in one step. Detrending the whole series first and deseasonalising afterwards does not do the same: a seasonal cycle much larger than the trend is not orthogonal to a linear fit over a finite window, and the leakage leaves a spurious sawtooth in the residual.

Missing values are not skipped: a NaN sample makes the fit of its calendar day NaN for that series, so all of that series' samples on that day become NaN.

source

AtmosTransport.Analysis.drift_statistics Method

drift_statistics(series) -> (; drift, increment_mean, increment_std, t, n)

Change of a global-total series $s_1,\dots ,s_N$ (for example a <tracer>_total_mass series, $N\ge 2$) relative to its initial value $s_1>0$:

• drift $=(s_N-s_1)/s_1$, the net relative change;

• the $n=N-1$ relative increments $d_k=(s_{k+1}-s_k)/s_1$, with their mean increment_mean $\overline{d}$ and sample standard deviation increment_std ${\sigma }_d$ (normalised by $n-1$);

• t $=\overline{d}/({\sigma }_d/\sqrt{n})$, the mean increment in units of its standard error.

Without sources or sinks the total is invariant, so drift is accumulated arithmetic error, and t tells a bias (a mean per-step error, $|t|\gtrsim 3$) from unbiased noise ($|t|\lesssim 2$). With uniformly spaced frames the increments are per snapshot interval.

Degenerate cases: with $n=1$ the standard deviation and t are NaN; with ${\sigma }_d=0$ (every increment equal), t is 0 when the increments are 0, as for a Float64 total that stays constant to the last bit, and $\pm \mathrm{\infty }$ otherwise.

Replaces increment_stats and the drift column of float32_conservation_attribution.py. (That script returned $t=0$ for any ${\sigma }_d=0$ or NaN.)

source

AtmosTransport.Analysis.expected_storage_change Function

expected_storage_change(method::AbstractExpectedEmission, ctx::RunContext,
                        spec::TransportTracerSpec, t0, t1) -> Float64

Tracer storage in kg that the tracer of spec gains from its sources between t0 and t1 (seconds since ctx.reference_time), by method. A tracer without a surface flux gains 0.

source

AtmosTransport.Analysis.gc_path Method

gc_path(archive, t::DateTime) -> String

Path of the archive's file for time t: dir/YYYY/MM/<name> if that file exists, otherwise dir/<name> (which need not exist), with <name> = GEOSChem.CATRINE_inst.YYYYMMDD_HHMMz.nc4.

source

AtmosTransport.Analysis.gc_times Method

gc_times(archive, start::DateTime, stop::DateTime) -> Vector{DateTime}

The archive's output times start, start + interval, … up to and including stop (when it falls on the grid).

source

AtmosTransport.Analysis.global_totals Method

global_totals(archive, times::AbstractVector{DateTime}) -> Vector

global_totals of each time, on all Julia threads (start Julia with -t N); an element is nothing where the file is unreadable. Each thread holds one file in memory while it is decoded.

source

AtmosTransport.Analysis.global_totals Method

global_totals(archive, t::DateTime) -> (; dry_air, burden, emission_rate) or nothing

Global sums of one output file of a GEOS-Chem archive, in GEOS-Chem's dry-air-equivalent storage units (kg of dry air carrying the tracer's mole fraction, the units of AtmosTransport's <tracer>_total_mass):

• dry_air $=\sum \mathrm{Met}\mathrm{\_}\mathrm{AD}$ [kg];

• burden[i] $=\sum {\chi }_i\mathrm{Met}\mathrm{\_}\mathrm{AD}$ [kg] for the dry mole fraction ${\chi }_i$ of tracer i of archive.tracers;

• emission_rate[i] $=(\sum E_i\mathrm{Met}\mathrm{\_}\mathrm{AREAM}2)M_{\mathrm{air}}/M_i$ [kg s⁻¹], the emission flux $E_i$ [kg species m⁻² s⁻¹] integrated over the sphere and converted with GEOS-Chem's molar masses ($M_{\mathrm{air}}$ = DRY_AIR_MOLAR_MASS, GEOS-Chem's AIRMW; $M_i$ the tracer's molar_mass). The emission fields are instantaneous: the rate at time t.

Products and sums are taken in Float64 over the Float32 fields; fill values count as zero. Returns nothing when the file is unreadable (recorded in the archive, see with_gc_file).

source

AtmosTransport.Analysis.logp_interpolate! Method

logp_interpolate!(out, p, p_mid, q; outside = NaN) -> out

In-place logp_interpolate for a vector of target pressures p.

source

AtmosTransport.Analysis.logp_interpolate Method

logp_interpolate(p, p_mid, q; outside = NaN)

Value of the column profile q at the pressure p, linear in log pressure.

p_mid holds the layer mid pressures of one column, strictly increasing (top first) and positive; q the layer values. p is one target pressure (the result is a Float64) or a vector of them (the result is a Vector{Float64}), in the unit of p_mid. For a target with $p_k\le p\le p_{k+1}$,

$$q(p)=q_k(1-w)+q_{k+1}w,w=\frac{\ln p-\ln p_k}{\ln p_{k+1}-\ln p_k}.$$

A target outside $[p_1,p_L]$ (above the top mid level or below the lowest one) gets outside: there is no extrapolation and no clamping, so a zonal mean on pressure levels averages only the columns that reach the level.

source

AtmosTransport.Analysis.nbands Method

nbands(bands::AbstractLatitudeBands) -> Int

Number of bands.

source

AtmosTransport.Analysis.report_dict Method

report_dict(report::BudgetReport) -> Dict{String, Any}

The JSON report of a budget check. The keys of check_trendy_npp_rh_run.py are kept (days, tracers, max_relative_budget_error, tolerance, passed, and per tracer under checks: expected_storage_change_kg, actual_storage_change_kg, relative_error, absolute_tolerance_kg (the strict tolerance) and passed); the other keys add the operational level, the carrier, the error as ppm of the carrier air, the column check and the skipped tracers. days is the budget interval in days.

source

AtmosTransport.Analysis.run_snapshot_paths Method

run_snapshot_paths(ctx::RunContext) -> Vector{String}

The snapshot files the run writes, in time order, named as the runner names them from [output].path: one file for split = "single"; for split = "daily" one file per transport binary that takes a snapshot (the binary's YYYYMMDD date fills {YYYYMMDD}, {date} or {day}). A snapshot at a window end belongs to that window's binary, the one at hour 0 to the first. Only NetCDF output can be read back; an ATMSNAP (format = "binary_mmap") or disabled output is an error.

source

AtmosTransport.Analysis.seasonal_trend_decomposition Method

seasonal_trend_decomposition(series, seen; min_months = 24,
                             min_years_per_month = 5, first_month = 1)
    -> (; anomaly, seasonal, detrended)

Split monthly series into a mean seasonal cycle, a linear trend and interannual anomalies, fitting only the observed months.

series has size (nloc, nmonth) or (nloc, nmonth, nseries): locations (latitude bands, regions or cells), consecutive calendar months starting in calendar month first_month, and optionally several series (observations and models) fitted on the same mask. seen (nloc × nmonth, Bool) marks the observed months; series must be finite wherever seen is true and is not read elsewhere.

For each location with at least min_months observed months, each series is reduced to its mean over the observed months, $v_k=x_k-\overline{x}$, and fitted on the observed months $k$ by

$$v_k\approx \sum _{m=1}^{12}{c}_m[\mathrm{month}(k)=m]+\beta t_k,t_k=\frac{k-(n+1)/2}{12}\text{years},$$

twelve calendar-month intercepts and a linear trend in time centred on the record, solved by the pseudoinverse (minimum norm, so a calendar month never observed gets $c_m=0$). Outputs, NaN wherever not computed:

• anomaly (size of series): the residual $v_k-{\hat{v}}_k$ on observed months;

• detrended (size of series): $v_k-\beta t_k$ on observed months, minus its mean over them;

• seasonal (nloc × 12 [× nseries], indexed by calendar month): $c_m$, set to NaN for calendar months observed in fewer than min_years_per_month years, then reduced by the mean of the remaining months.

Adding a constant to a series leaves all three outputs unchanged, and the decomposition of a sum of series is the sum of their decompositions.

source

AtmosTransport.Analysis.skill_metrics Method

skill_metrics(model, reference; min_n = 3, degenerate_std = 1e-12)
    -> (; r, rmse, bias, amplitude_ratio, n)

Unweighted agreement of model with reference (arrays of equal size) over the n pairs where both are finite:

• r: Pearson correlation; NaN when n < min_n or either population standard deviation is below degenerate_std (in the units of the data);

• rmse $=\sqrt{\frac{1}{n}\sum (x_i-y_i{)}^2}$;

• bias $=\frac{1}{n}\sum (x_i-y_i)$, model minus reference;

• amplitude_ratio $={\sigma }_x/{\sigma }_y$, population standard deviations (${\sigma }^2=\frac{1}{n}\sum (x_i-\overline{x}{)}^2$).

With n = 0 every statistic is NaN.

source

AtmosTransport.Analysis.twin_relative_difference Method

twin_relative_difference(candidate, reference) -> (; raw, net)

Frame-by-frame relative difference of a global total from its twin: candidate and reference hold the same total (for example <tracer>_total_mass) at the same frames, from two runs that differ in one setting, typically a Float32 run and its Float64 twin.

• raw[k] $=(c_k-r_k)/r_k$, NaN where $r_k=0$;

• net[k] $=$ raw[k] - raw[1] when $r_1\ne 0$: the difference accumulated after the first frame. The first frame carries the rounding of the initial state to the run precision (Float32(400e-6) is 2.53e-8 below 400e-6), a constant offset that is removed. When $r_1=0$, for a tracer that starts from zero mass, net = raw, relative to the reference total, which is then the mass emitted since the first frame.

Both series must start at the same frame, the start of the compared span: net is net of whatever difference that frame holds, so for the offset of the initial state it must be the initial snapshot of both runs. Non-finite totals are an error.

Replaces f32_minus_f64 of float32_conservation_attribution.py (which is net at the last common frame; its $(c_k-r_k)/(r_k-r_1)$ for a zero start is raw[k] exactly) and the rel/net series of compare_float32_float64_runs.py, with the same Float64 operations.

source

AtmosTransport.Analysis.unreadable_files Method

unreadable_files(archive) -> Vector{Pair{DateTime, String}}

The times whose files with_gc_file could not read so far, in time order, each with the reason (the path is gc_path).

source

AtmosTransport.Analysis.weighted_difference_stats Method

weighted_difference_stats(model, reference, weights) -> (; mean, rms, n)

Weighted mean and root mean square of the difference $d_i=x_i-y_i$ of model ($x$) and reference ($y$), for example a Float32 column-mean field and its Float64 twin with cell areas as weights. Only the n cells where both values are finite enter; their weights must be finite and non-negative and are normalised, $w_i\leftarrow w_i/\sum w$:

• mean $=\sum w_i{d}_i$;

• rms $=\sqrt{\sum w_i{d}_i^2}$.

The difference is formed per cell before any sum, and the sums are compensated (Neumaier) in linear cell order, so a mean far below the spread of the differences (a small bias under noise) keeps its relative accuracy. (The bias of weighted_fit_stats is a difference of two weighted means and loses the digits the two fields share.) When no cell has finite values, or their weights sum to 0, mean and rms are NaN.

Replaces the area-weighted daily mean and RMS of the column-mean difference in compare_float32_float64_runs.py, which let NaN propagate instead of leaving those cells out; on fields without NaN the two agree to the rounding of NumPy's pairwise sum (measured ≤ 1.5e-14 relative on C90 fields).

source

AtmosTransport.Analysis.weighted_fit_stats Method

weighted_fit_stats(model, reference, weights)
    -> (; r2, slope, bias, rmse, ref_mean, ref_std)

Weighted agreement of model ($y$) with reference ($x$), for example two fields on one grid with cell areas as weights. Only the cells where both values are finite enter; their weights must be finite and non-negative and are normalised, $w_i\leftarrow w_i/\sum w$. With the weighted means $\overline{x}=\sum wx$, $\overline{y}=\sum wy$ and moments $s_{xx}=\sum w(x-\overline{x}{)}^2$, $s_{yy}$, $s_{xy}$:

• r2 $=s_{xy}^2/(s_{xx}{s}_{yy})$, NaN unless both moments are positive;

• slope $=s_{xy}/s_{xx}$, the regression of the model on the reference (NaN when $s_{xx}=0$);

• bias $=\overline{y}-\overline{x}$;

• rmse $=\sqrt{\sum w(y-x{)}^2}$;

• ref_mean $=\overline{x}$ and ref_std $=\sqrt{s_{xx}}$.

When no cell has finite values (or the selected weights sum to 0) every statistic is NaN.

source

AtmosTransport.Analysis.with_gc_file Method

with_gc_file(f, archive, t::DateTime)

Open the archive's file for time t (gc_path) and return f(ds), where ds is the NCDataset; return nothing if the file cannot be read.

The file is read into memory first (without the lock, so threads read files concurrently), then opened from memory and passed to f under the package's NetCDF lock (Output.with_netcdf_lock; netcdf-c is not thread-safe). f should only copy the arrays it needs out of ds; computing on them after with_gc_file returns keeps the lock short. Memory use is the file size plus whatever f returns.

The file is unreadable when it is missing or cannot be read (SystemError, EOFError), when it is not a NetCDF/HDF5 file or is truncated (a NetCDFError on opening), and when f hits a NetCDFError (a missing variable, a failed read) or a variable without a time record: then the time and the reason are recorded in the archive (see unreadable_files) and nothing is returned. Any other error propagates, including a KeyError or BoundsError raised by f itself.

source