Post Processing Functions

gcmprocpy provides a range of functions for post processing the data. All derived variables are registered in a central registry and dispatched automatically when passed to any plot function — no special handling is required.

Note

For live examples with output, see the Emissions notebook.

Derived Variables

gcmprocpy supports derived variables — quantities computed from multiple dataset fields rather than read directly. These are used like any other variable name in plot and data extraction functions.

Category	Variable Name	Description	Requirements
Emissions	`NO53`	5.3-micron NO emission	Temperature, Atomic Oxygen, NO
Emissions	`CO215`	15-micron CO2 emission	Temperature, Atomic Oxygen, CO2
Emissions	`OH83`	OH v(8,3) band emission (simple)	Temperature, Atomic Oxygen, O2, N2
OH Meinel	`OH_<upper>_<lower>`	Specific OH Meinel band (e.g. `OH_8_3`)	Temperature, O, O2, N2, H, O3, HO2
OH Meinel	`OH_TOTAL`	Sum of all 39 Meinel band emissions	Temperature, O, O2, N2, H, O3, HO2
OH Meinel	`OH_VIB_<v>`	Vibrational level population (v=0-9)	Temperature, O, O2, N2, H, O3, HO2
EP Flux	`EPVY`	Meridional EP flux component	Temperature, U, V winds
EP Flux	`EPVZ`	Vertical EP flux component	Temperature, U, V, W winds
EP Flux	`EPVDIV`	EP flux divergence (wave forcing)	Temperature, U, V, W winds

Example: Using derived variables in plot functions

datasets = gy.load_datasets(directory, dataset_filter)
time = '2022-01-01T12:00:00'

# Emissions — use with plt_lat_lon
plot = gy.plt_lat_lon(datasets, 'NO53', time=time, level=4.0)
plot = gy.plt_lat_lon(datasets, 'CO215', time=time, level=4.0)
plot = gy.plt_lat_lon(datasets, 'OH_8_3', time=time, level=4.0)

# EP flux — use with plt_lev_lat
plot = gy.plt_lev_lat(datasets, 'EPVY', time=time)
plot = gy.plt_lev_lat(datasets, 'EPVDIV', time=time)

# Species density — use arr_density for explicit unit control
cm3 = gy.arr_density(datasets, 'O1', time=time, to_unit='CM3')

Emissions

Simple Emissions Plots

gcmprocpy provides functions for computing airglow and infrared emissions from model output.

Note

The emission rate equations operate on number densities (cm⁻³). Species inputs are converted automatically via the density machinery, reading each field’s units attribute — so emissions are correct whether a history stores species as mass mixing ratio (kg/kg), volume mixing ratio (mol/mol), or number density (cm-3). N2 is derived as 1 − O2 − O when absent. (TIE-GCM thermosphere histories typically lack the mesospheric chemistry — H, O3, HO2 — the full OH model needs.)

Example 1: Plotting 5.3 micron NO emission

datasets = gy.load_datasets(directory, dataset_filter)
variable_name = 'NO53'
time = '2022-01-01T12:00:00'
pressure_level = 4.0
plot = gy.plt_lat_lon(datasets, variable_name, time=time, level=pressure_level)

Example 2: Plotting 15 micron CO2 emission

datasets = gy.load_datasets(directory, dataset_filter)
variable_name = 'CO215'
time = '2022-01-01T12:00:00'
pressure_level = 4.0
plot = gy.plt_lat_lon(datasets, variable_name, time=time, level=pressure_level)

Example 3: Plotting OH emission for the v(8,3) band

datasets = gy.load_datasets(directory, dataset_filter)
variable_name = 'OH83'
time = '2022-01-01T12:00:00'
pressure_level = 4.0
plot = gy.plt_lat_lon(datasets, variable_name, time=time, level=pressure_level)

Emissions Array Functions

5.3 micron NO emission

This function processes the given datasets to generate an array of 5.3-micron NO emissions based on temperature, O1, and NO data.

gcmprocpy.data_emissions.arr_mkeno53(datasets, variable_name, time, selected_lev_ilev=None, selected_unit=None, plot_mode=False)[source]

Generate 5.3-micron NO emission array from datasets. This function processes the given datasets to generate an array of 5.3-micron NO emissions based on temperature, O1, and NO data.

Requires the following variables: - TN or T: Temperature - O1 or O: Oxygen - NO: Nitric oxide concentration

Parameters:

datasets (list) – List of datasets to process.
variable_name (str) – Name of the variable to process.
time (datetime) – Specific time for which data is to be processed.
selected_lev_ilev (int, optional) – Selected level or ilev. Defaults to None.
selected_unit (str, optional) – Unit of the variable. Defaults to None.
plot_mode (bool, optional) – Flag to indicate if plot mode is enabled. Defaults to False.

Returns:

If plot_mode is False, returns an numpy array containing 5.3-micron NO emissions for the specified time and level. If plot_mode is True, returns a tuple containing: - NO_emission (ndarray): Array of 5.3-micron NO emissions. - level (ndarray): Array of levels. - unique_lats (ndarray): Array of unique latitudes. - unique_lons (ndarray): Array of unique longitudes. - str: Empty string placeholder. - str: Description of the emission (“5.3-micron NO”). - selected_mtime (datetime): Selected time. - filename (str): Name of the file processed.

Return type:

Union[xarray.DataArray, tuple]

15 micron CO2 emission

This function processes the given datasets to generate an array of 15-micron CO2 emissions based on temperature, O1, and CO2 data.

gcmprocpy.data_emissions.arr_mkeco215(datasets, variable_name, time, selected_lev_ilev=None, selected_unit=None, plot_mode=False)[source]

Generate CO2 emissions using the mkeco215 model based on temperature, oxygen, and CO2 data.

Requires the following variables: - TN or T: Temperature - O1 or O: Oxygen - CO2: CO2 concentration

Parameters:

datasets (list) – List of datasets to be used for extracting variables.
variable_name (str) – Name of the variable to be processed.
time (datetime) – Specific time for which the data is to be processed.
selected_lev_ilev (int, optional) – Specific level or ilevel to be selected. Defaults to None.
selected_unit (str, optional) – Unit to which the data should be converted. Defaults to None.
plot_mode (bool, optional) – If True, additional plotting information is returned. Defaults to False.

Returns:

If plot_mode is False, returns an numpy array containing CO2 emissions for the specified time and level. If plot_mode is True, returns a tuple containing: - CO2_emission (numpy.ndarray): CO2 emissions calculated by the mkeco215 model. - level (int): Selected level or ilevel. - unique_lats (numpy.ndarray): Unique latitudes. - unique_lons (numpy.ndarray): Unique longitudes. - str: Placeholder string. - str: Long name for the 15-micron CO2 emission. - datetime: Selected time. - str: Filename of the dataset.

Return type:

Union[xarray.DataArray, tuple]

OH emission for the v(8,3) band

This function processes the given datasets to generate an array of OH emissions for the v(8,3) band based on temperature, O1, and OH data.

gcmprocpy.data_emissions.arr_mkeoh83(datasets, variable_name, time, selected_lev_ilev=None, selected_unit=None, plot_mode=False)[source]

Generate OH emissions using the mkeoh83 model based on temperature, oxygen, and nitrogen data.

Requires the following variables: - TN or T: Temperature - O1 or O: Oxygen - O2: Molecular oxygen - N2: Molecular nitrogen

Parameters:

datasets (list) – List of datasets to be used for extracting variables.
variable_name (str) – Name of the variable to be processed.
time (datetime) – Specific time for which the data is to be processed.
selected_lev_ilev (int, optional) – Specific level or ilevel to be selected. Defaults to None.
selected_unit (str, optional) – Unit to which the data should be converted. Defaults to None.
plot_mode (bool, optional) – If True, additional plotting information is returned. Defaults to False.

Returns:

If plot_mode is False, returns a numpy array containing OH emissions for the specified time and level. If plot_mode is True, returns a tuple containing: - OH_emission (numpy.ndarray): OH emissions calculated by the mkeoh83 model. - level (int): Selected level or ilevel. - unique_lats (numpy.ndarray): Unique latitudes. - unique_lons (numpy.ndarray): Unique longitudes. - str: Placeholder string. - str: Long name for the OH v(8,3) emission. - datetime: Selected time. - str: Filename of the dataset.

Return type:

Union[numpy.ndarray, tuple]

Raw Emissions Calculations

These are the underlying physics functions that compute emission rates from raw arrays. They can be called directly for custom processing pipelines.

5.3 micron NO emission

Calculates 5.3 micron NO emission (from John Wise).

gcmprocpy.data_emissions.mkeno53(arr_temp, arr_o, arr_no)[source]

Calucates 5.3 micron NO emission (from John Wise).

Ported from tgcmproc mkemiss.F (subroutine mkeno53). cgs units: Tk in K, [O]/[NO] as number density (cm⁻³), result in photons cm⁻³ s⁻¹. (The Fortran used a truncated pi=3.14156; here np.pi.)

The formula used is:

N(5.3 mic) = (2.63E-22) * exp[-2715 / Tk] * [O] * [NO]
             -----------------------------------------
               (4 * Pi) * (10.78 + 6.5E-11 * [O])

Where:

[O] is the oxygen concentration.
[NO] is the nitric oxide concentration.
Tk is the temperature in Kelvin.

\[N(5.3 \, \mu m) = \frac{2.63 \times 10^{-22} \cdot \exp\left(-\frac{2715}{T_k}\right) \cdot [O] \cdot [NO]}{4 \pi \cdot \left(10.78 + 6.5 \times 10^{-11} \cdot [O]\right)}\]

Parameters:

arr_temp (numpy.ndarray) – Array of temperatures in Kelvin.
arr_o (numpy.ndarray) – Array of oxygen concentrations.
arr_no (numpy.ndarray) – Array of nitric oxide concentrations.

Returns:

Calculated NO emission at 5.3 microns.

Return type:

numpy.ndarray

15 micron CO2 emission

Calculates 15 micron CO2 emission (from John Wise).

gcmprocpy.data_emissions.mkeco215(arr_temp, arr_o, arr_co2)[source]

Calucates 15 micron CO2 emission (from John Wise).

Ported from tgcmproc mkemiss.F (subroutine mkeco215). cgs units: Tk in K, [O]/[CO2] as number density (cm⁻³), result in photons cm⁻³ s⁻¹. (The Fortran used a truncated pi=3.14156; here np.pi.)

The formula used is:

N(15 mic) = (5.94E-26) * sqrt(Tk) * exp[-960 / Tk] * [O] * [CO2]
             ---------------------------------------------------
                 (4 * Pi) * (1.28 + 3.5E-13 * sqrt(Tk) * [O])

Where:

[O] is the oxygen concentration.
[CO2] is the carbon dioxide concentration.
Tk is the temperature in Kelvin.

The 15 micron term is only the O-CO2 collisional component, but it accounts for at least 80% of the radiance above 110 km.

\[N(15 \, \mu m) = \frac{5.94 \times 10^{-26} \cdot \sqrt{T_k} \cdot \exp\left(-\frac{960}{T_k}\right) \cdot [O] \cdot [CO_2]}{4 \pi \cdot \left(1.28 + 3.5 \times 10^{-13} \cdot \sqrt{T_k} \cdot [O]\right)}\]

Parameters:

arr_temp (numpy.ndarray) – Array of temperatures (Tk).
arr_o (numpy.ndarray) – Array of oxygen concentrations [O].
arr_co2 (numpy.ndarray) – Array of CO2 concentrations [CO2].

Returns:

Calculated 15 micron CO2 emission.

Return type:

numpy.ndarray

OH emission for the v(8,3) band

Calculates OH emission for the v(8,3) band.

gcmprocpy.data_emissions.mkeoh83(arr_temp, arr_o, arr_o2, arr_n2)[source]

Calculate OH emission for the v(8,3) band.

Ported from tgcmproc mkemiss.F (subroutine mkeoh83, “6/95: Return Oh emission v(8,3) band”):

fout = f8*[O]*[O2]*(pk6n2*[N2] + pk6o2*[O2]) / (260 + 2e-11*[O2])

with f8=0.29, pk6n2=5.70e-34*(300/T)**2.62, pk6o2=5.96e-34*(300/T)**2.37. cgs units: T in K, densities in cm⁻³, emission in photons cm⁻³ s⁻¹.

Parameters:

arr_temp (numpy.ndarray) – Array of temperatures (K).
arr_o (numpy.ndarray) – Array of atomic oxygen densities (cm^-3).
arr_o2 (numpy.ndarray) – Array of molecular oxygen densities (cm^-3).
arr_n2 (numpy.ndarray) – Array of molecular nitrogen densities (cm^-3).

Returns:

Calculated OH emission for the v(8,3) band.

Return type:

numpy.ndarray

OH Meinel Band Model

The full OH Meinel band vibrational emission model solves coupled steady-state rate equations for 10 vibrational levels (v=0 through v=9) and computes emission rates for all 39 Meinel bands. Ported from tgcmproc ohrad.F (subroutine ohrad; B. Foster, U. B. Makhlouf, SDL/Stewart Radiance Lab).

Variable names use the pattern OH_<upper>_<lower> for specific bands (e.g. OH_8_3), OH_VIB_<v> for vibrational populations, and OH_TOTAL for the total emission rate.

Note

The full OH model requires seven species in the dataset: temperature, O, O2, N2, H, O3, and HO2.

Example 1: Plot OH(8,3) band emission

datasets = gy.load_datasets(directory, dataset_filter)
plot = gy.plt_lat_lon(datasets, 'OH_8_3', time='2022-01-01T12:00:00', level=4.0)

Example 2: Plot total OH emission

datasets = gy.load_datasets(directory, dataset_filter)
plot = gy.plt_lat_lon(datasets, 'OH_TOTAL', time='2022-01-01T12:00:00', level=4.0)

Example 3: Extract OH data for custom analysis

datasets = gy.load_datasets(directory, dataset_filter)

# Get PlotData for a specific band
result = gy.arr_mkoh_band(datasets, 'OH_8_3', time='2022-01-01T12:00:00',
                          selected_lev_ilev=4.0, plot_mode=True)
print(result.values.shape, result.variable_unit)

# Get raw array
values = gy.arr_mkoh_band(datasets, 'OH_8_3', time='2022-01-01T12:00:00',
                          selected_lev_ilev=4.0)

OH Array Function

gcmprocpy.data_oh.arr_mkoh_band(datasets, variable_name, time, selected_lev_ilev=None, selected_unit=None, plot_mode=False)[source]

Compute OH vibrational emission from datasets using the full model.

The variable_name determines what is returned:

'OH_<upper>_<lower>' — a specific band, e.g. 'OH_8_3'
'OH_VIB_<v>' — vibrational level population, e.g. 'OH_VIB_8'
'OH_TOTAL' — sum of all 39 band emission rates

Requires these variables in the datasets: temperature (TN/T), O2, O (O1), N2, H, O3, HO2.

Parameters:

datasets (list) – List of ModelDataset objects.
variable_name (str) – Encoded output selector (see above).
time – Time value to extract.
selected_lev_ilev – Pressure level or 'mean'.
selected_unit (str, optional) – Desired unit (unused, kept for API compat).
plot_mode (bool) – If True, return a PlotData object.

Returns:

numpy.ndarray or PlotData.

Raises:

ValueError – If required species are missing or variable_name is invalid.

OH Physics Calculation

The underlying steady-state solver for vibrational populations and band emission rates.

gcmprocpy.data_oh.ohrad(temp, o2, o, n2, h, o3, ho2, oh=None)[source]

Full OH Meinel band vibrational emission model.

Port of tgcmproc ohrad.F (B. Foster, U. B. Makhlouf, SDL). Solves the coupled steady-state rate equations for 10 vibrational levels (v = 0 … 9) at each grid point independently.

The steady-state system at each point is A x = b where A is a 10×10 rate matrix (radiative decay, collisional quenching, multi-quantum O2 redistribution, reverse quenching) and b is the external production vector (H + O3 primary, O + HO2 secondary).

References

Einstein A coefficients — Nelson et al. (1990), updated Makhlouf (1999) Quenching rates — Adler-Golden (1997) Production branching — Klenerman & Smith (H+O3), Kaye (O+HO2)

The 10×10 rate matrix A and production vector b mirror tgcmproc ohrad.F (subroutine ohrad; Fortran DECOMP/SOLVE → numpy.linalg.solve()); band rate = [OH(v)] · A(v→v') in photons cm⁻³ s⁻¹ (the Fortran iwatts=0 branch). With oh=None the v=0 self-recycling term 2·rk8·[OH] is omitted (tgcmproc always supplies it from the model).

Parameters:

temp – Temperature (K). Any shape (...).
o2 – O2 number density (cm⁻³), same shape.
o – Atomic oxygen number density (cm⁻³), same shape.
n2 – N2 number density (cm⁻³), same shape.
h – Atomic hydrogen number density (cm⁻³), same shape.
o3 – Ozone number density (cm⁻³), same shape.
ho2 – HO2 number density (cm⁻³), same shape.
oh – OH ground-state density (cm⁻³), same shape. Optional — if None, the v = 0 recycling terms are omitted.

Returns:

vib_pop (ndarray, shape (..., 10)): Vibrational populations [OH(v)] in cm⁻³.
band_emission (dict): {(upper_v, lower_v): emission_array} in photons cm⁻³ s⁻¹, each array of shape (...).

Return type:

tuple

Eliassen-Palm Flux

Eliassen-Palm (EP) flux diagnostics quantify wave-mean flow interaction in the atmosphere. gcmprocpy computes three components from 3-D wind and temperature fields, ported from tgcmproc epflux.F (subroutines epfluxy/epfluxz/epfluxdiv; B. Foster and Hanli Liu, 1998):

EPVY — meridional EP flux component (m² s^-2)
EPVZ — vertical EP flux component (m² s^-2)
EPVDIV — EP flux divergence / wave forcing (m s^-1 day^-1)

EP flux variables produce level-latitude cross sections and are plotted with plt_lev_lat. EPVY requires only horizontal winds; EPVZ and EPVDIV additionally require vertical wind (W/OMEGA).

Example 1: Plot EP flux divergence

datasets = gy.load_datasets(directory, dataset_filter)
plot = gy.plt_lev_lat(datasets, 'EPVDIV', time='2022-01-01T12:00:00')

Example 2: Plot meridional EP flux

datasets = gy.load_datasets(directory, dataset_filter)
plot = gy.plt_lev_lat(datasets, 'EPVY', time='2022-01-01T12:00:00')

Example 3: Extract EP flux data for custom analysis

datasets = gy.load_datasets(directory, dataset_filter)

result = gy.arr_epflux(datasets, 'EPVDIV', time='2022-01-01T12:00:00')
print(result.values.shape)   # (nlev, nlat)
print(result.variable_unit)  # 'm s⁻¹ day⁻¹'

EP Flux Array Function

gcmprocpy.data_epflux.arr_epflux(datasets, component, time, log_level=True)[source]

Compute an EP flux component from model datasets.

Extracts temperature and wind fields at the requested time, converts units to SI, and calls epflux().

Parameters:

datasets (list) – List of ModelDataset objects.
component (str) – 'EPVY', 'EPVZ', or 'EPVDIV'.
time – Timestamp (numpy.datetime64 or ISO string).
log_level (bool) – Whether to log-transform level coordinates for the output PlotData.

Returns:

Result with shape (nlev, nlat) suitable for plt_lev_lat().

Return type:

PlotData

Note

EPVDIV needs the zonal-mean total mass density rhozm. Following tgcmproc mkrhokg (mkderived.F), it is built by summing the major species (O, O₂, N₂) converted to mass density via arr_density() (to_unit='GM/CM3', ×1000 → kg m⁻³), so it carries the pkt = p/(k_B·T) vertical falloff — for both TIE-GCM and WACCM-X and for any source unit (MMR / VMR / cm⁻³). If the required species (or, for WACCM-X, the hybrid-pressure coefficients hyam/hybm/PS) are unavailable, it falls back to the ideal-gas proxy ρ ∝ exp(-lev)/T.

Raises:: ValueError – If component is invalid or required variables are missing from the datasets.

EP Flux Physics Calculation

The core computation that takes 3-D wind/temperature arrays and returns zonal-mean EP flux components. Can be called directly with numpy arrays for custom workflows.

gcmprocpy.data_epflux.epflux(temp, u, v, lats, levs, w=None, rho=None)[source]

Compute Eliassen-Palm flux components.

All inputs must be in SI units (K, m s⁻¹, degrees, dimensionless log-pressure levels).

The calculation follows the quasi-geostrophic EP flux formulation:

\[ \begin{align}\begin{aligned}S_y = -\overline{u'v'} + \frac{\overline{v'T'}}{\gamma} \frac{\partial \bar{u}}{\partial z}\\S_z = -\left[\overline{u'w'} + \left(\frac{1}{\cos\phi} \frac{\partial(\bar{u}\cos\phi)}{\partial y} - f\right) \frac{\overline{v'T'}}{\gamma}\right]\\D = \frac{1}{\cos^2\phi} \frac{\partial(\cos^2\phi \, S_y)}{\partial y} + \frac{1}{\bar{\rho}} \frac{\partial(\bar{\rho} \, S_z)}{\partial z}\end{aligned}\end{align} \]

where primes denote deviations from the zonal mean, overbars are zonal means, γ is the static stability, and f the Coriolis parameter.

Ported from tgcmproc epflux.F (module epflux; subroutines epfluxy/epfluxz/epfluxdiv with zonal-mean, gamma and zpht setup in save_epv; B. Foster & H. Liu, 2-4/98). EPVY, EPVZ, γ, the reference scale height H = r·ts/g, and the 3/12/98 sign fix (Sy = -u'v' + (v'T'/γ)·∂u/∂z, plus sign) follow the Fortran. The EPVDIV mass density is supplied by arr_epflux() (see its note).

Parameters:

temp – Temperature (K), shape (nlev, nlat, nlon).
u – Zonal wind (m s⁻¹), same shape.
v – Meridional wind (m s⁻¹), same shape.
lats – Latitude array (degrees), shape (nlat,).
levs – Log-pressure levels (dimensionless), shape (nlev,).
w – Vertical wind (m s⁻¹), same shape as temp. Optional — required for EPVZ and EPVDIV.
rho – Total mass density (kg m⁻³), same shape as temp. Optional — used as zonal mean for EPVDIV. If omitted, falls back to the ideal-gas proxy ρ ∝ exp(-lev)/T (constant mean molecular mass), which deviates from tgcmproc when M̄ varies with altitude (notably in the thermosphere).

Returns:

Keys 'EPVY', 'EPVZ', 'EPVDIV', each a (nlev, nlat) ndarray. EPVZ and EPVDIV are None if w is not provided.

Return type:

dict

Species-Aware Density Conversions

gcmprocpy converts atmospheric species fields among the four density representations used by tgcmproc, ported from denconv.F (subroutines mkdenparms/denconv; B. Foster):

MMR — mass mixing ratio (dimensionless, mass fraction)
CM3 — number density (molecules cm^-3)
CM3-MR — volume mixing ratio / mole fraction (dimensionless)
GM/CM3 — mass density (g cm^-3)

Cross-unit conversions require the mean air molar mass (BARM) and the air number density (pkt), both derived from the model’s temperature and O / O₂ fields. Formulas:

BARM = 1 / (O₂/32 + O/16 + (1 − O₂ − O)/28)
pkt  = pressure / (k_B · T)                    [CGS]

MMR → CM3     :   f · pkt · BARM / W
MMR → CM3-MR  :   f · BARM / W
MMR → GM/CM3  :   f · pkt · BARM · 1.66e-24

The source unit is read from each field’s units attribute (kg/kg → MMR, mol/mol → CM3-MR, cm-3 → CM3, …), so conversion works whether a history stores species as mass mixing ratio, volume mixing ratio, or number density.

Note

Both TIE-GCM and WACCM-X are supported. The only model-specific step is how pressure is recovered from the vertical coordinate: TIE-GCM uses log-pressure p = p₀·exp(−ζ); WACCM-X uses the CAM hybrid sigma-pressure p = hyam·P0 + hybm·PS (so the file must carry hyam/hybm/PS).

Example: Convert atomic oxygen from MMR to number density

import gcmprocpy as gy

datasets = gy.load_datasets(directory, dataset_filter)

# Reads O1 from the dataset, uses the 'units' attr for the source unit
# (falls back to passing from_unit=... explicitly)
result = gy.arr_density(datasets, 'O1', time='2022-01-01T12:00:00',
                         to_unit='CM3')
print(result.values.shape)      # (nlev, nlat, nlon)
print(result.variable_unit)     # 'CM3'

Example: Direct conversion with explicit barm / pkt

from gcmprocpy import (
    convert_density_units, compute_barm, compute_pkt,
    get_species_molar_mass,
)

# Using raw numpy arrays — barm/pkt can be scalar or full 3-D fields
barm = compute_barm(o_mmr=o1_field, o2_mmr=o2_field)
pkt  = compute_pkt(levs, temperature, model='TIE-GCM')
w    = get_species_molar_mass('TIE-GCM', 'O2')

o2_cm3 = convert_density_units(o2_mmr, 'MMR', 'CM3',
                                barm=barm, pkt=pkt, molar_mass=w)

Density Conversion Array Function

gcmprocpy.data_density.arr_density(datasets, variable_name, time, *, to_unit='CM3', from_unit=None, selected_lev_ilev=None, log_level=True, **kwargs)

Core Conversion and Helpers

The pure-math functions that arr_density() composes. Use these directly for custom workflows (e.g. converting in-memory arrays from external sources).

gcmprocpy.data_density.convert_density_units(values, from_unit, to_unit, *, barm, pkt, molar_mass)[source]

Convert a density field between MMR / CM3 / CM3-MR / GM/CM3.

Uses MMR as the pivot. Ported from tgcmproc denconv.F (subroutine denconv, module den_convert):

MMR→CM3    = mmr * pkt * barm / w
MMR→CM3-MR = mmr * barm / w
MMR→GM/CM3 = mmr * pkt * barm * 1.66e-24
CM3→MMR    = f   * w / (pkt * barm)

where barm is the mean air molar mass and pkt the air number density (compute_barm() / compute_pkt()), and w the species molar mass. The reverse (non-MMR → MMR) directions are algebraic inversions of the above (not written verbatim in the Fortran).

Parameters:

values – Array-like field to convert.
from_unit – Source / target unit strings. Accepts common aliases ('cm-3', 'kg/kg', 'g/cm3'…).
to_unit – Source / target unit strings. Accepts common aliases ('cm-3', 'kg/kg', 'g/cm3'…).
barm – Mean air molar mass (g/mol) broadcast-compatible with values.
pkt – Air number density (cm⁻³), same broadcast shape.
molar_mass – Species molar mass (g/mol), scalar.

Returns:

ndarray with the same shape as values, in to_unit.

gcmprocpy.data_density.compute_barm(o_mmr, o2_mmr)[source]

Mean air molar mass (g/mol) from atomic-O and molecular-O₂ MMRs.

Ported from tgcmproc denconv.F (subroutine mkdenparms, module den_convert, standard/TIE-GCM branch):

barm = 1 / (o2/32. + o1/16. + max(.00001, 1.-o2-o1) / 28.)

The residual (1 − O₂ − O) is treated as N₂ (molar mass 28); the .00001 floor is the original Fortran guard for the upper thermosphere where O₂ + O ≈ 1. (The planetary jtgcm/mtgcm/vtgcm barm variants are not ported.)

gcmprocpy.data_density.compute_pkt(levs, temperature, model='TIE-GCM', *, hyam=None, hybm=None, p0=None, ps=None)[source]

Air number density (cm⁻³) = pressure / (k_B · T), in CGS.

Ported from tgcmproc denconv.F (subroutine mkdenparms; pkt formula at line 68):

pkt(:,k) = p0 * exp(-(zpb + (k-1)*dz)) / (boltz * TN),   boltz = 1.3805e-16

The TIE-GCM log-pressure coordinate gives p = p0·exp(-ζ) with p0 = 5e-4 (cgs; proc.F labels it “mb” but the value matches the microbar = dyn cm⁻² convention). The WACCM-X hybrid-pressure branch below is a gcmprocpy extension with no tgcmproc counterpart.

The pressure is recovered from the model’s vertical coordinate, which differs by model:

TIE-GCM (log-pressure lev = ζ = ln(p₀/p)):
```
p = p₀ · exp(−ζ),   p₀ = 5 × 10⁻⁴ dyn cm⁻²
```
temperature lev axis must be axis 0 when multi-dimensional; lev broadcasts across the remaining axes.
WACCM-X (CAM hybrid sigma-pressure):
```
p(k,j,i) = hyam(k) · P0 + hybm(k) · PS(j,i)     [Pa]
```
converted to dyn cm⁻² (×10). Requires hyam/hybm (on the field’s vertical coordinate) and surface pressure ps (lat, lon); p0 defaults to 1 × 10⁵ Pa if not supplied. Returns a full (nlev, nlat, nlon) field since pressure varies horizontally.

Parameters:

levs – Vertical coordinate values, shape (nlev,).
temperature – Temperature in K.
model – 'TIE-GCM' or 'WACCM-X'.
hyam – Hybrid A/B coefficients on the field’s level coordinate (WACCM-X only), shape (nlev,).
hybm – Hybrid A/B coefficients on the field’s level coordinate (WACCM-X only), shape (nlev,).
p0 – Reference pressure in Pa (WACCM-X only; default 1e5).
ps – Surface pressure in Pa, shape (nlat, nlon) (WACCM-X only).

gcmprocpy.data_density.get_species_molar_mass(model, variable_name)[source]

Molar mass (g/mol) for a species variable name in the given model.

Resolves through get_species_names() so that model-specific naming (e.g. TIE-GCM O1 vs WACCM-X O) is handled uniformly.

Raises:: ValueError – If the variable is not a known species in the model, or no molar mass is registered for that species role.

Supported unit aliases include 'cm-3' → CM3, 'kg/kg' → MMR, 'mol/mol' → CM3-MR, and 'g/cm3' → GM/CM3. See SUPPORTED_DENSITY_UNITS for the canonical tuple.

Difference Fields

gcmprocpy supports computing difference fields between two datasets (e.g. perturbed vs control runs). Both raw differences and percent differences are supported.

Example: Compute and plot a difference field

datasets_pert = gy.load_datasets(pert_dir, dataset_filter)
datasets_ctrl = gy.load_datasets(ctrl_dir, dataset_filter)
time = '2022-01-01T12:00:00'

pert_result = gy.arr_lat_lon(datasets_pert, 'TN', time=time,
                             selected_lev_ilev=4.0, plot_mode=True)
ctrl_result = gy.arr_lat_lon(datasets_ctrl, 'TN', time=time,
                             selected_lev_ilev=4.0, plot_mode=True)

# Raw difference (perturbation - control)
diff_result = gy.diff_plotdata(pert_result, ctrl_result, diff_type='RAW')

# Percent difference
diff_result = gy.diff_plotdata(pert_result, ctrl_result, diff_type='PERCENT')

gcmprocpy.data_diff.compute_diff(pert_values, cntr_values, diff_type='RAW')[source]

Compute difference between perturbation and control arrays.

Parameters:

pert_values (numpy.ndarray) – Perturbation field values.
cntr_values (numpy.ndarray) – Control field values.
diff_type (str) – ‘RAW’ for simple subtraction, ‘PERCENT’ for percentage difference (pert - cntr) / cntr * 100.

Returns:

Difference values. For PERCENT mode, grid points where |cntr| <= 1e-20 are set to zero to avoid division by zero.

Return type:

numpy.ndarray

Raises:

ValueError – If pert_values and cntr_values have different shapes or diff_type is not ‘RAW’ or ‘PERCENT’.

gcmprocpy.data_diff.diff_plotdata(pert_result, cntr_result, diff_type='RAW')[source]

Compute difference between two PlotData objects.

Parameters:

pert_result (PlotData) – Perturbation result from an arr_* function.
cntr_result (PlotData) – Control result from the same arr_* function.
diff_type (str) – ‘RAW’ or ‘PERCENT’.

Returns:

New PlotData with difference values and an updated variable_long_name that includes the diff type label.

Return type:

PlotData

Raises:

ValueError – If the two PlotData value arrays have different shapes.

Derived Variable Registry

All derived variables (emissions, OH bands, EP flux) are managed through a central registry. When a derived variable name is passed to any plot function, the registry automatically dispatches to the correct handler. You can also register custom derived variables.

gcmprocpy.containers.get_species_names(model)[source]

Return species name mapping for a model type.

Uses MODEL_DEFAULTS as the single source of truth for mapping canonical role names to dataset variable names.

Parameters:: model (str) – Model type ('TIE-GCM' or 'WACCM-X').
Returns:: Mapping from canonical names (e.g. 'temp', 'o', 'o2') to dataset variable names (e.g. 'TN', 'O1', 'O2').
Return type:: dict
Raises:: ValueError – If model is not recognized.

gcmprocpy.containers.register_derived(name, handler, plot_types=None)[source]

Register a derived variable computation handler.

Parameters:

name (str) – Variable name (e.g. 'NO53') or a glob-style pattern ending with * (e.g. 'OH_*').
handler (callable) – Function with signature (datasets, variable_name, time, **kwargs) -> PlotData.
plot_types (set, optional) – Plot types this variable supports (e.g. {'lat_lon', 'lev_lat'}). None means all.

gcmprocpy.containers.resolve_derived(variable_name)[source]

Look up the handler for a derived variable name.

Checks exact matches first, then pattern matches (keys ending with *).

Parameters:: variable_name (str) – The variable name to look up.
Returns:: (handler, True) if found, (None, False) otherwise.
Return type:: tuple

Example: Check if a variable is derived

from gcmprocpy import resolve_derived

handler, is_derived = resolve_derived('EPVY')
print(is_derived)   # True

handler, is_derived = resolve_derived('TN')
print(is_derived)   # False

Example: Register a custom derived variable

from gcmprocpy import register_derived

def my_custom_var(datasets, variable_name, time, **kwargs):
    # Custom computation ...
    return result

register_derived('MY_VAR', my_custom_var)

Derivable Fields (auto-derive if missing)

In addition to the derived output registry above, gcmprocpy has a lower-level derivable-field layer for intermediate quantities the model often does not write to the history file. When a plot or data-extraction function asks for a field that is absent from the dataset, it is computed on the full native grid from the fields that are present and injected transparently — so it can be sliced and plotted like any stored variable. A real field in the file always takes priority over the derived form.

Currently registered derivables:

Name	Definition	Inputs
`N2`	molecular nitrogen mixing-ratio residual `max(1e-5, 1 - O2 - O)`	O2, O
`O/N2` (alias `O_N2`)	atomic oxygen / molecular nitrogen ratio	O, N2
`N2/O` (alias `N2_O`)	molecular nitrogen / atomic oxygen ratio	N2, O
`O/O2` (alias `O_O2`)	atomic oxygen / molecular oxygen ratio	O, O2
`O/O2+N2` (alias `O_O2pN2`)	atomic oxygen / (O2 + N2) ratio	O, O2, N2
`RHO`	total air mass density `O2 + O + N2` (g cm⁻³)	TN, O, O2, N2
`PMB`	pressure in millibars (from the model vertical coordinate)	TN (+ `PS` for WACCM-X)
`TNFP`	frost-point temperature (K) — requires `H2O`	TN, O, O2, H2O
`OX`	odd oxygen `O + O3` (additive group, members’ unit)	O, O3
`NOZ`	odd nitrogen `NO + NO2`	NO, NO2
`HOX`	odd hydrogen `OH + HO2 + H`	OH, HO2, H
`O/CO2` (alias `O_CO2`)	atomic oxygen / CO2 ratio	O, CO2

This is what lets the full OH model and the OH83 emission run on histories that omit N2 (it is derived from O2/O automatically). Derivables can also depend on other derivables (e.g. O/N2 builds on the derived N2). The RHO/PMB/TNFP fields reuse the density machinery, so they are unit-aware (MMR / VMR / cm⁻³) and work for both TIE-GCM and WACCM-X; TNFP raises a chain-aware error on histories without H2O. The composition groups OX/NOZ/HOX are summed from their members in the file’s shared unit (exact — densities are additive; tgcmproc writes the group field natively and gcmprocpy reconstructs it when absent), require the members to share a unit, and raise a chain-aware error when a member is absent (e.g. OH/NO2 are not on standard TIE-GCM histories).

Register a custom derivable:

from gcmprocpy import register_derivable

# formula receives a dict of input DataArrays keyed by canonical role
register_derivable('O/O', lambda inp, mds: inp['o'] / inp['o2'],
                    inputs=['o', 'o2'], units='ratio',
                    long_name='O / O2 ratio')

gcmprocpy.containers.register_derivable(name, formula, inputs, units='', long_name=None, models=None)[source]

Register a derivable intermediate field (computed when absent from a file).

Parameters:

name (str) – Output field name (the name a user/handler requests, e.g. 'N2' or 'O/N2').
formula (callable) – fn(inputs_by_role: dict, mds) -> xarray.DataArray operating on full-grid DataArrays.
inputs (list[str]) – Canonical role names ('o2', 'o', 'n2', …) or literal field names this field is built from.
units (str) – Units to stamp on the result.
long_name (str) – Long name; defaults to name.
models (set, optional) – Restrict to these model types; None = any.

gcmprocpy.containers.resolve_derivable(variable_name)[source]: Return the DERIVABLE_VARIABLES entry for variable_name, or None.

Persisting Derived Quantities to NetCDF

Calculated derivable fields can be written back into their source NetCDF file so that subsequent loads read them directly instead of recomputing. The variable is computed on the full grid and appended in place to the file.

import gcmprocpy as gy

datasets = gy.load_datasets('/path/to/output', dataset_filter='sech')
gy.save_derived(datasets, ['N2', 'O/N2'])   # append once
# next session: N2 and O/N2 are already in the file, read like any variable

Note

This appends to the original file in place, so the file must be writable — read-only archive histories raise PermissionError (persist into a writable copy instead). NetCDF cannot delete a variable in place, so a field already present is skipped. Only derivable intermediate fields are persisted this way; slice-based derived outputs (emissions, OH bands, EP flux) are not.

gcmprocpy.io.save_derived(datasets, variable_names, overwrite=False, verbose=True)[source]

Compute derivable field(s) on the full native grid and append them in place to each dataset’s source NetCDF file, so subsequent loads read them directly instead of recomputing.

Only derivable intermediate fields are persisted this way — quantities computed on the full grid from other fields (e.g. 'N2' and the composition ratios). Slice-based derived outputs (emissions, OH bands, EP flux) are not handled here.

The dataset is closed, the variable is appended via netCDF4, and the file is reopened, so the in-memory datasets keep working and the data cache is cleared.

Parameters:

datasets (list[ModelDataset]) – Loaded datasets (modified in place).
variable_names (str | list[str]) – Field name(s) to compute and write.
overwrite (bool) – NetCDF cannot delete a variable in place, so a field already present on disk is skipped with a warning regardless; regenerate into a fresh copy to replace it.
verbose (bool) – Log progress.

Returns:

"<path>:<var>" entries actually written.

Return type:

list[str]

Raises:

ValueError – If a name is not present and not derivable, or a dataset has no on-disk source path.
PermissionError – If a source file is read-only or locked.

tgcmproc Provenance & References

The physics routines documented above are ports of NCAR’s legacy tgcmproc Fortran post-processor. Each formula was cross-referenced against its original source and verified for numerical equivalence; the table records the provenance — source file, subroutine, and original author/attribution — so every ported formula can be traced back to the Fortran.

gcmprocpy function	tgcmproc source	Attribution	Notes
`mkeno53` (`NO53`)	`mkemiss.F` :: `mkeno53`	John Wise	5.3 µm NO emission; cgs number densities → photons cm⁻³ s⁻¹
`mkeco215` (`CO215`)	`mkemiss.F` :: `mkeco215`	John Wise	15 µm CO₂ emission (O–CO₂ collisional term)
`mkeoh83` (`OH83`)	`mkemiss.F` :: `mkeoh83`	tgcmproc (6/95)	OH v(8,3) band emission
`ohrad` (OH Meinel model)	`ohrad.F` :: `ohrad`	Foster, U. B. Makhlouf (SDL/Stewart Radiance Lab)	10-level steady-state OH(v) model; 39 Meinel bands
`epflux` (`EPVY`/`EPVZ`/`EPVDIV`)	`epflux.F` :: `epfluxy`/`epfluxz`/`epfluxdiv`/`save_epv`	Foster & Hanli Liu (1998)	Eliassen–Palm flux & divergence; see EPVDIV caveat below
`compute_barm`	`denconv.F` :: `mkdenparms`	Foster	mean air molar mass `1/(o2/32+o1/16+max(.00001,1-o2-o1)/28)` (TIE-GCM branch)
`compute_pkt`	`denconv.F` :: `mkdenparms`	Foster	air number density `p0·exp(-ζ)/(boltz·TN)`; WACCM-X hybrid-pressure branch is a gcmprocpy extension
`convert_density_units`	`denconv.F` :: `denconv`	Foster	MMR ↔ CM3 / CM3-MR / GM-CM3 (reverse paths are algebraic inversions)
`get_species_molar_mass`	`fset_known.F` :: `fset_known`	tgcmproc	species molar masses (`flds_known(n)%wt`)
`N2` derivable	`mkderived.F` :: `mkderived`	tgcmproc	residual `fn2 = max(.00001, 1-o2-o1)`
`O/N2`, `N2/O`, `O/O2`, `O/O2+N2` ratios	`mkderived.F` :: `mkderived`	tgcmproc	composition ratios; unit-invariant, formed on the source fields
`RHO` derivable	`mkderived.F` :: `mkderived` (L222) / `mkrhokg`	tgcmproc	total mass density `O2+O+N2` → GM/CM3 (carries the `pkt` falloff)
`PMB` derivable	`mkderived.F` :: `mkderived` (L798)	tgcmproc	`(n_O2+n_O+n_N2)·k_B·T·1e-3` = `p·1e-3` (mb)
`TNFP` derivable	`mkderived.F` :: `mkderived` (L803)	tgcmproc	frost point `T − 6077.4/(28.548 − ln(p_H2O[mb]))`; requires `H2O`
`OX` / `NOZ` / `HOX` groups	`fset_known.F` (native fields; `fcomponents` + group `%wt` 16/30/17)	tgcmproc	odd O/N/H groups (O+O3, NO+NO2, OH+HO2+H); summed from members when absent
`O/CO2` ratio	`mkderived.F` :: `mkderived` (L618)	tgcmproc	atomic oxygen / CO2 ratio (terrestrial branch); unit-invariant

Scientific references

EP flux — Volland, H. (1988), Atmospheric Tidal and Planetary Waves, Kluwer Academic Publishers, §2.7.
OH Meinel model — Einstein A coefficients: Nelson et al. (1990), updated Makhlouf (1999); quenching rates: Adler-Golden (1997); production branching: Klenerman & Smith (H + O₃), Kaye (O + HO₂).
NO 5.3 µm / CO₂ 15 µm emissions — formulations from John Wise (via tgcmproc mkemiss.F).

Note

EPVDIV mass density. Following tgcmproc mkrhokg (mkderived.F), arr_epflux builds the zonal-mean mass density for EPVDIV by summing the major species (O, O₂, N₂) converted to GM/CM3 via arr_density (×1000 → kg m⁻³), so it carries the pkt = p/(k_B·T) vertical falloff for both TIE-GCM and WACCM-X and for any source unit (MMR / VMR / cm⁻³). If the required species — or, for WACCM-X, the hybrid-pressure coefficients (hyam/hybm/PS) — are unavailable, it falls back to the ideal-gas proxy ρ ∝ exp(-lev)/T.