API Reference

Every public function and class in pvtend is documented below with full parameter descriptions. Click [source] next to any entry to jump directly to the implementation.

Grid & Constants

The grid module provides the standard 1.5° Northern-Hemisphere grid used throughout the package, spatial cropping/interpolation utilities, and event-centred patch extraction. Constants include Earth parameters, thermodynamic constants, and default level lists.

`NHGrid`(lat, lon)	Northern Hemisphere regular lat-lon grid.
`default_nh_grid`()	Return the standard 1.5° NH grid (90°N–0°, -180°–180°).
`EventPatch`(grid[, lat_half, lon_half])	Event-centred patch extraction from a full NH grid.
`constants`	Physical and numerical constants for pvtend.

Derivatives

Centred finite-difference operators on lat-lon-pressure grids. All operators handle:

Periodic zonal boundaries — the globe wraps in longitude (toggled via periodic flag).
One-sided differences at polar / top / bottom boundaries.
Non-uniform pressure spacing for ddp().

ddx, ddy accept both 2-D (nlat, nlon) and 3-D (nlev, nlat, nlon) inputs; ddp operates along axis 0; ddt uses centred differences with one-sided at the temporal boundaries.

`ddx`(field, dx_arr[, periodic])	Zonal derivative ∂f/∂x with periodic wrapping.
`ddy`(field, dy)	Meridional derivative ∂f/∂y with one-sided differences at boundaries.
`ddp`(field, plevs_pa)	Pressure derivative ∂f/∂p for non-uniform levels.
`ddt`(field, dt_s)	Time derivative ∂f/∂t via centred differences.
`d_dlambda`(field, dlambda[, periodic])	Zonal derivative in angular coordinates ∂f/∂λ.
`d_dphi`(field, dphi)	Meridional derivative in angular coordinates ∂f/∂φ.
`div_spherical`(Ax, Ay, lat_rad, dphi, dlambda)	Spherical divergence of a 2-D vector field (Aλ, Aφ).

Spherical-harmonic operators

Since v2.10, pvtend ships a spherical-harmonic backend (pvtend.sh_ops) built on pyspharm that is the new default for the Helmholtz decomposition, the QG ω velocity-potential inversion, and mixed second derivatives. All public operators accept Northern- Hemisphere arrays (nlat_nh, nlon) (lat = 0…90°) and internally mirror them to the global sphere via parity_mirror() (scalars and u use even reflection, v uses odd) before the SH projection — this avoids the polar singularities of pure finite differences and yields machine-precision round-trip Helmholtz identities on band-limited inputs.

`gradient_sh`(field, lat, lon[, R_earth])	Physical gradient of a scalar via spherical harmonics.
`laplacian_sh`(field, lat, lon[, R_earth])	Spherical Laplacian of a scalar field via spherical harmonics.
`invert_laplacian_sh`(rhs, lat, lon[, ...])	Solve `∇²χ = rhs` on the sphere with the global-mean zero gauge.
`vortdiv_sh`(u, v, lat, lon[, R_earth])	Vorticity ζ and divergence δ from a 2-D wind field.
`helmholtz_sh`(u, v, lat, lon[, R_earth, ...])	Spherical-harmonic Helmholtz decomposition.
`second_derivs_sh`(field, lat, lon[, R_earth])	Second-order horizontal derivatives `(f_xx, f_xy, f_yy)`.
`parity_mirror`(field, lat, kind)	Mirror an NH-only field onto the full globe with the correct parity.
`restrict_to_nh`(field_global, lat_global)	Slice a full-globe field to the Northern-Hemisphere portion.
`is_nh_grid`(lat[, tol])	Detect whether lat covers only the Northern Hemisphere.

Climatology

Compute and load an hourly climatology for ERA5 pressure-level variables. The climatology is smoothed in time (day-of-year) using a low-pass Fourier filter (4 harmonics) and optionally with 2 zonal modes, producing per-variable, per-month NetCDF files for anomaly computation.

The Helmholtz climatology pre-computes the rotational/divergent decomposition of the climatological wind for each month (24 NetCDF files). These are loaded during the tendency pipeline to compute anomaly Helmholtz fields by subtraction.

`compute_climatology`(data_dir, output_dir, years)	Compute per-variable per-month hourly climatology.
`load_climatology`(clim_path[, engine, ...])	Load climatology, auto-detecting file layout.
`compute_helmholtz_climatology`(clim_dir, ...)	Pre-compute Helmholtz decomposition of the climatological wind.
`load_helmholtz_climatology`(clim_dir, month, *)	Load pre-computed Helmholtz-decomposed climatological wind fields.

QG Omega Equation

Solves the Hoskins (1978) quasi-geostrophic omega equation:

\[\nabla^2_p \omega + \frac{f_0^2}{\sigma}\,\frac{\partial^2\omega} {\partial p^2} \;=\; -2\,\nabla\cdot\mathbf{Q} - \beta\,\frac{R_d}{\sigma\,p}\,\frac{\partial T}{\partial x}\]

where \(\mathbf{Q}\) is the Q-vector (Hoskins et al. 1978).

Two solver methods

log20 (default) — Strongly Implicit Procedure (SIP, Stone 1968) with a full 3-D spherical finite-difference stencil in conservative cos(φ:sub:`j±½`) flux-divergence form (Lynch 1989; cf. xinvert). This replaces the previous tan φ coefficient form (v ≤ 2.10) which was ill-conditioned at high latitudes; the polar metric now vanishes naturally as cos(φ_j±½) → 0. Closest analogue to Li & O’Gorman (2020). Numba-accelerated (nogil=True) when available (~3–6 s per event); falls back to pure-Python if numba is not installed. Always solved on the full Northern Hemisphere grid (periodic in longitude, Neumann at equatorial and polar faces) and the event patch is extracted afterward. At φ = ±90° the polar Dirichlet row is replaced by the zonal-mean of omega_b so that only the geometrically well-defined \(m=0\) mode survives at the geometric pole — matching the Helmholtz spherical-FFT and SH solvers, where high-\(m\) modes vanish through the \(\cos^{-2}\varphi\) eigenvalue scaling.
sp19 — Steinfeld & Pfahl (2019) empirical scaling. Sets \(\omega_\text{dry} = \frac{1}{3}\,\omega_\text{total}\). Zero-cost, preserves spatial structure. QG-moist ω equals the moist residual (no separate term-C solve).

Both methods apply a latitude taper (linearly 0 → 1 between 15°N and 25°N, tapering back to 0 above 85°N) to all RHS terms (A, B, and C) to suppress polar metric-term singularities and enforce QG validity.

Boundary conditions (LOG20)

Inhomogeneous Dirichlet: real ERA5 \(\omega\) is prescribed on all six faces (top, bottom, and four lateral boundaries) via the omega_b parameter. This anchors the QG inversion to the reanalysis vertical velocity field at the domain edges, eliminating the bias of homogeneous (\(\omega = 0\)) boundaries.

Diabatic heating — term C

The solver supports two independent diabatic forcings, and is called three times per timestep:

LOG20 term C (budget-closure, Li & O’Gorman 2020): When the Eulerian temperature tendency \(\partial T/\partial t\) is available (dT_dt parameter), the diabatic heating rate is diagnosed from the thermodynamic equation:

\[J = \underbrace{c_p\!\left(\frac{\partial T}{\partial t} + \mathbf{v}\cdot\nabla_p T\right)}_{J_1} + \underbrace{\left(-\frac{\sigma\,p}{R_d}\;c_p\;\omega\right)}_{J_2}\]

\[\text{Term C}_\text{LOG20} = -\frac{\kappa}{p}\;\nabla^2_s\,J\]

This captures all diabatic processes (LHR + radiation + sensible heat flux + diffusion).

Emanuel term \(C_\text{em}\) (parameterised LHR, Emanuel 1987; Tamarin & Kaspi 2016): The latent heat release is parameterised from thermodynamic profiles during saturated ascent (\(\omega < 0\)):

\[\dot\theta_\text{LHR} = \omega\!\left(\frac{\partial\theta}{\partial p} - \frac{\gamma_m}{\gamma_d}\;\frac{\theta}{\theta_E}\; \frac{\partial\theta_E}{\partial p}\right)\]

\[J_\text{em} = c_p\;\dot\theta_\text{LHR}\;\frac{T}{\theta}, \qquad C_\text{em} = -\frac{\kappa}{p}\;\nabla^2_s\,J_\text{em}\]

This captures only latent heating during saturated ascent.

Three-solve strategy:

QG-dry (terms A+B, homogeneous top/bottom BCs) → \(\omega_\text{dry}\)
QG-total (A+B+\(C_\text{LOG20}\), ERA5 BCs) → \(\omega_\text{qg,total}\)
Emanuel (A+B+\(C_\text{em}\), ERA5 BCs) → \(\omega_\text{em}\)

From these, four omega components are derived:

\[\begin{split}\omega_\text{qg\_moist} &= \omega_\text{qg,total} - \omega_\text{dry} \\ \omega_\text{em\_moist} &= \omega_\text{em} - \omega_\text{dry} \\ \omega_\text{moist} &= \omega_\text{total} - \omega_\text{dry}\end{split}\]

This yields a 4-way vertical-velocity decomposition:

\[\omega = \underbrace{\omega_\text{dry}}_{\text{A+B}} + \underbrace{\omega_\text{qg\_moist}}_{\text{all diabatic (LOG20)}} + \underbrace{(\omega_\text{moist} - \omega_\text{qg\_moist})}_{\text{non-QG residual}}\]

with the Emanuel-moist component \(\omega_\text{em\_moist}\) providing a sub-partition that isolates the LHR-only contribution. Comparing \(\omega_\text{qg\_moist}\) and \(\omega_\text{em\_moist}\) separates LHR from non-LHR diabatic effects.

References: Hoskins B J, Draghici I, Davies H C (1978) Q.J.R.M.S. 104, 31–38. Li L, O’Gorman P A (2020) J. Climate. Emanuel K A, Fantini M, Thorpe A J (1987) J. Atmos. Sci. 44, 1559–1573. Tamarin T, Kaspi Y (2016) J. Atmos. Sci. 73, 1687–1707. Stone H L (1968) SIAM J. Numer. Anal. 5, 530–558. Steinfeld D, Pfahl S (2019) Clim. Dyn. 53, 6159–6180.

solve_qg_omega_sip(u, v, t, lat, lon, ...[, ...])

Solve the QG omega equation using the SIP (Strongly Implicit Procedure).

ω blowup detection

ERA5 occasionally exhibits unphysically large vertical velocities at upper levels (typically thunderstorm artefacts in the 4D-Var increment) that contaminate downstream tendency budgets. pvtend flags such timestamps with a fixed Pa s⁻¹ hard cutoff at the chosen pressure level — the canonical rule is \(|\omega| > 5\,\text{Pa s}^{-1}\) at 300 hPa.

Earlier versions of pvtend silently clipped ω at ±5 Pa/s inside solve_qg_omega_sip(), which corrupted the affected events; that behaviour is OFF by default in v2.10+. Detection is now performed post-hoc by scan_omega_blowup(), which emits a CSV of offending timestamps for downstream exclusion by the event-NPZ builders, RWB classifier, and composite stage.

scan_omega_blowup(era5_w_glob[, level_pa, ...])

Flag ERA5 timestamps whose max absolute ω at level_pa exceeds threshold.

CLI:

pvtend blowup-scan \
    --era5 '/net/flood/data2/users/x_yan/era/era5_w_*.nc' \
    --level 300 --threshold 5.0 \
    --out outputs/blowup_scan/omega_300hPa_5pa.csv

Helmholtz Decomposition

Decomposes a 2-D wind field \((u, v)\) into three orthogonal components:

\[\mathbf{u} = \underbrace{\mathbf{u}_\text{rot}}_{\text{non-divergent}} + \underbrace{\mathbf{u}_\text{div}}_{\text{irrotational}} + \underbrace{\mathbf{u}_\text{har}}_{\text{harmonic}}\]

where the rotational wind derives from a streamfunction \(\psi\) (\(u_\text{rot}=-\partial\psi/\partial y\), \(v_\text{rot}=\partial\psi/\partial x\)), the divergent wind from a velocity potential \(\chi\) (\(u_\text{div}=\partial\chi/\partial x\), \(v_\text{div}=\partial\chi/\partial y\)), and the harmonic residual satisfies \(\nabla\cdot\mathbf{u}_\text{har}=0\) and \(\zeta(\mathbf{u}_\text{har})=0\).

The decomposition proceeds by:

Computing spherical vorticity \(\zeta\) and divergence \(\delta\) from the input wind via compute_vorticity_divergence():

\[\zeta = \frac{1}{a\cos\varphi}\!\left[ \frac{\partial v}{\partial\lambda} - \frac{\partial(u\cos\varphi)}{\partial\varphi}\right] = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} + \frac{u\tan\varphi}{a}\]

\[\delta = \frac{1}{a\cos\varphi}\!\left[ \frac{\partial u}{\partial\lambda} + \frac{\partial(v\cos\varphi)}{\partial\varphi}\right] = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} - \frac{v\tan\varphi}{a}\]

The \(\pm\tan\varphi/a\) metric terms arise from spherical curvature and are significant at blocking latitudes (30–70°N).
Removing the area-weighted (\(\cos\varphi\)) mean from each (Fredholm solvability condition).
Solving Poisson equations \(\nabla^2\psi=\zeta\) and \(\nabla^2\chi=\delta\) using the spherical FFT Poisson solver solve_poisson_spherical_fft().
Recovering \((u_\text{rot},v_\text{rot})\) and \((u_\text{div},v_\text{div})\) via the spectral gradient gradient().
Computing the harmonic residual by subtraction.

Spherical Poisson solver

The Poisson equations are solved using a conservative spherical Laplacian (following MiniUFO/xinvert) that combines:

FFT in longitude (exploiting 360° periodicity), and
a tridiagonal solve in latitude with the full \(\cos\varphi\) metric factors.

The solver supports both Dirichlet (default, for backward compatibility) and Neumann (ghost-point) boundary conditions in latitude via the bc_type parameter. Since v2.1, helmholtz_decomposition() uses bc_type="neumann" with physical BCs derived from the ERA5 wind field, reducing the harmonic residual from ~44% to ~5% on the bounded NH domain.

References: Lynch P (1989) MWR 117, 1492–1500.

Spectral gradient (wind recovery)

After solving the Poisson equations, the divergent and rotational winds are recovered via gradient(). The zonal derivative uses a spectral (FFT) method:

\[\frac{\partial\chi}{\partial x}\bigg|_{\!j} = \frac{1}{\Delta x_j}\;\mathrm{IFFT}\!\left( \mathrm{i}\,\frac{2\pi m}{N}\,\widehat{\chi}_m\right)\]

where \(\widehat{\chi}_m\) are the Fourier coefficients along a latitude circle of \(N\) points. This ensures that the discrete composition \(\nabla\cdot(\nabla\chi)\) is consistent with the compact spherical Laplacian used by the Poisson solver, eliminating the \(\cos^2(m\pi/N)\) attenuation that centred finite differences would introduce. The Nyquist mode is zeroed for even \(N\).

The meridional derivative retains centred finite differences (the \(\varphi\) direction is non-periodic).

`helmholtz_decomposition`(u, v, lat, lon[, ...])	Helmholtz decomposition on a lat-lon grid.
`helmholtz_decomposition_3d`(u_3d, v_3d, lat, lon)	Helmholtz decomposition for each vertical level.

Solver backends

Both helmholtz_decomposition() and helmholtz_decomposition_3d() accept a method keyword that selects the underlying numerical backend:

`method`	Description
`"spectral"`	Default — spherical-harmonic transform via `pvtend.sh_ops` (pyspharm). NH arrays are parity-mirrored to the global sphere before the SH projection. Machine-precision round-trip identities, no FD truncation error, no polar metric residual.
`"sh"`	Synonym for `"spectral"`.
`"fft"`	Spherical-FFT Poisson solver (`pvtend.helmholtz._helmholtz_decomposition_fd`) used in pvtend ≤ 2.9; combines an FFT in longitude with a tridiagonal solve in latitude. Retained for legacy reproducibility.
`"fd"`	Synonym for `"fft"`.

The method="spectral" (SH) path is the recommended default — it produces a harmonic residual ≲ 1 % of \(\|\mathbf{u}\|\) on band-limited NH inputs and reduces the previous 5 % residual of the FFT path.

Low-level Poisson helpers (importable from pvtend.helmholtz):

`solve_poisson_spherical_fft`(rhs, lat, dy, ...)	Solve the spherical Poisson equation — periodic in λ, configurable φ BCs.
`laplacian_spherical_fft`(phi, lat, dy, dlon_rad)	Apply the spherical Laplacian using the SAME stencil as the Poisson solver.
`compute_vorticity_divergence`(u, v, dx, dy[, ...])	Spherical vorticity and divergence.
`gradient`(phi, dx, dy)	Physical gradient (∂φ/∂x, ∂φ/∂y) — spectral zonal, FD meridional.

laplacian_spherical_fft is the forward operator matching the Poisson solver’s conservative spherical Laplacian stencil. It enables machine-precision round-trip verification: \(\mathcal{L}(\text{solve\_poisson}(f)) = f\) (interior points).

Moist/Dry Omega

Partitions total vertical velocity into four components and recovers the corresponding divergent winds:

\[\omega = \underbrace{\omega_\text{dry}}_{\text{QG (A+B)}} + \underbrace{\omega_\text{qg\_moist}}_{\text{QG (A+B+C}_\text{LOG20}\text{) − (A+B)}} + \underbrace{\omega_\text{residual}}_{\text{non-QG}}\]

with \(\omega_\text{moist} = \omega_\text{total} - \omega_\text{dry}\) encompassing both QG-moist and non-QG contributions.

An additional Emanuel-moist component \(\omega_\text{em\_moist} = \omega_{A+B+C_\text{em}} - \omega_{A+B}\) isolates the LHR-only response (see QG Omega Equation for details).

The divergent-wind recovery proceeds as follows:

\(\omega_\text{moist} = \omega_\text{total} - \omega_\text{dry}\)
Solve \(\nabla^2\chi_\text{moist} = -\partial\omega_\text{moist}/\partial p\) on each pressure level (spherical Poisson solver).
\((u_\text{div,moist}, v_\text{div,moist}) = \nabla\chi_\text{moist}\)
Solve \(\nabla^2\chi_\text{dry} = -\partial\omega_\text{dry}/\partial p\) independently on each pressure level (same spherical Poisson solver).
\((u_\text{div,dry}, v_\text{div,dry}) = \nabla\chi_\text{dry}\)
QG-moist divergent wind via the same Poisson inversion applied to \(\omega_\text{qg\_moist}\): \((u_\text{div,qg\_moist}, v_\text{div,qg\_moist}) = \nabla\chi_\text{qg\_moist}\)
Emanuel-moist divergent wind via the same Poisson inversion applied to \(\omega_\text{em\_moist}\): \((u_\text{div,em\_moist}, v_\text{div,em\_moist}) = \nabla\chi_\text{em\_moist}\)

All four divergent-wind components are recovered via independent Poisson inversions using solve_chi_from_omega(). This avoids amplification of discretisation errors from the Helmholtz decomposition that would occur with a residual-based approach (\(\mathbf{u}_{\text{div,dry}} = \mathbf{u}_{\text{div}} - \mathbf{u}_{\text{div,moist}}\)).

Equation	Solver	Function
\(\nabla^2\psi = \zeta\) (Helmholtz)	Spherical FFT Poisson	`solve_poisson_spherical_fft()`
\(\nabla^2\chi = \delta\) (Helmholtz)	Spherical FFT Poisson	`solve_poisson_spherical_fft()`
\(\nabla^2\chi_\text{moist} = -\partial\omega_\text{moist}/\partial p\)	Spherical FFT Poisson	`solve_chi_from_omega()`
\(\nabla^2\chi_\text{dry} = -\partial\omega_\text{dry}/\partial p\)	Spherical FFT Poisson	`solve_chi_from_omega()`
\(\nabla^2\chi_\text{qg\_moist} = -\partial\omega_\text{qg\_moist}/\partial p\)	Spherical FFT Poisson	`solve_chi_from_omega()`
\(\nabla^2\chi_\text{em\_moist} = -\partial\omega_\text{em\_moist}/\partial p\)	Spherical FFT Poisson	`solve_chi_from_omega()`
\(\mathcal{L}(\chi) = f\) (forward Laplacian for verification)	Conservative spherical stencil	`laplacian_spherical_fft()`

Note: the horizontal wind decomposition remains 2-way (dry + moist). The QG-moist and Emanuel-moist divergent winds are physically interpretable sub-partitions of the full moist divergent wind — they can substitute for \(\mathbf{u}_{\chi,\text{moist}}\) in PV cross-term analyses but are not additional additive components of the Helmholtz decomposition. Comparing their PV cross-terms isolates LHR from non-LHR diabatic effects.

Total-field approximation for horizontal divergence. The QG omega solve and Poisson inversion are performed on total fields (\(\omega\), not \(\omega'\)), exploiting the dominance of the anomaly vertical velocity in midlatitude synoptic systems (\(|\omega'| \gg |\bar\omega|\)). Because the Poisson operator is linear, the same approximation propagates to the horizontal divergent wind:

\[\omega_\text{moist} \approx \omega'_\text{moist} \;\;\Longrightarrow\;\; \mathbf{u}_{\chi,\text{moist}} \approx \mathbf{u}'_{\chi,\text{moist}}\]

This avoids a costly second anomaly-field QG inversion while capturing the dominant moist–dry partition of both vertical velocity and horizontal divergent wind.

`decompose_omega`(omega_total, ...[, u_div, v_div])	Full adiabatic/diabatic omega decomposition with divergent wind recovery.
`solve_chi_from_omega`(omega, lat, lon, plevs_pa)	Solve ∇²χ = -∂ω/∂p at each level.
`verify_div_additivity`(u_div, ...)	Verify additive consistency of divergent wind decomposition.

Orthogonal Basis Decomposition

Projects PV tendency fields onto a set of six orthogonal basis functions constructed via Gram-Schmidt orthogonalisation. Each basis captures a distinct dynamical mode of blocking intensification / weakening:

\(\Phi_1 = q'\) — PV anomaly (intensification / weakening)
\(\Phi_2 = \partial q/\partial x\) — Zonal gradient (zonal propagation)
\(\Phi_3 = \partial q/\partial y\) — Meridional gradient (meridional propagation)
\(\Phi_4 = \partial^2 q/\partial x\,\partial y\) — Shear deformation (\(\gamma_1\))
\(\Phi_5 = \partial^2 q/\partial x^2 - \partial^2 q/\partial y^2\) — Normal-strain deformation (\(\gamma_2\))
\(\Phi_6 = \partial^2 q/\partial x^2 + \partial^2 q/\partial y^2\) — Laplacian dispersion (\(\sigma\))

The projection coefficients \(\beta, a_x, a_y, \gamma_1, \gamma_2, \sigma\) quantify each term’s contribution to the blocking lifecycle at each timestep.

Because the raw basis fields span many orders of magnitude (PVU, PVU m^-1, PVU m^-2), each is scaled to \(O(1)\) before orthogonalisation. The default constants are:

Pre-normalization constants
Basis	Constant	Typical magnitude
\(\Phi_1\) (PV anomaly)	\(10^{6}\)	\(\sim 10^{-6}\) K m² s^-1 kg^-1
\(\Phi_2, \Phi_3\) (gradients)	\(10^{12}\)	\(\sim 10^{-12}\) K m s^-1 kg^-1
\(\Phi_4, \Phi_5, \Phi_6\) (2nd derivatives)	\(10^{18}\)	\(\sim 10^{-18}\) K s^-1 kg^-1

Alternatively, auto_prenorm() computes per-event normalization as \(1/\mathrm{median}(|f|)\) over the finite masked region, activated via prenorm_mode="auto".

When the mask threshold produces multiple disconnected regions (common at large Eulerian displacement), compute_orthogonal_basis() automatically retains only the single connected component that encloses the patch centre. If the centre falls outside all blobs, the blob whose nearest boundary pixel is closest to the centre is selected. This uses scipy.ndimage.label for connected-component analysis.

compute_orthogonal_basis() supports temporal down-scaling from hourly ERA5 snapshots to a sub-hourly evaluation instant via built-in bi-linear interpolation. When the pv_anom_prev, pv_dx_prev, pv_dy_prev, and (optionally) pv_dx_dy_prev, pv_dx_dx_prev, pv_dy_dy_prev keyword arguments are supplied (the fields at hour dh − 1), the positional fields (at hour dh) are interpolated before basis construction:

\[f_{\alpha} = \alpha\,f_{dh} + (1 - \alpha)\,f_{dh-1}\]

The default \(\alpha = 1.0\) uses only the current-time fields (no interpolation). Set \(\alpha = 0.75\) to place the evaluation instant at dh − 15 min (¾ of the way from dh − 1 toward dh), effectively down-scaling from 1-hour to 15-minute temporal resolution.

Because linear interpolation commutes with spatial differentiation on a fixed grid, the interpolated gradients satisfy \(\nabla f_{\alpha} = \alpha\,\nabla f_{dh} + (1-\alpha)\,\nabla f_{dh-1}\), and the cross-derivative \(\partial^2 q/\partial x\,\partial y\) inherits the same property.

Temporal interpolation parameters
Parameter	Description
`pv_anom_prev`	PV anomaly \(q'\) at dh − 1 (the earlier snapshot).
`pv_dx_prev`	Zonal PV gradient \(\partial q/\partial x\) at dh − 1.
`pv_dy_prev`	Meridional PV gradient \(\partial q/\partial y\) at dh − 1.
`pv_dx_dy_prev`	Cross derivative \(\partial^2 q/\partial x\,\partial y\) at dh − 1.
`pv_dx_dx_prev`	\(\partial^2 q/\partial x^2\) at dh − 1.
`pv_dy_dy_prev`	\(\partial^2 q/\partial y^2\) at dh − 1.
`interp_alpha`	Weight for the positional (current-dh) fields (default 1.0). Set to 0.75 for 15-min before dh, or 0 to use only dh − 1.

All three _prev fields must be provided together (or none). When omitted the function falls back to the original single-snapshot behaviour with no interpolation.

A standalone helper is also available for manual interpolation outside the basis pipeline:

lerp_fields(fields_prev, fields_curr[, ...])

Linearly interpolate selected 2-D fields between two time-steps.

`OrthogonalBasisFields`(phi_int, phi_dx, ...)	Container for the six orthogonal basis fields.
`compute_orthogonal_basis`(pv_anom, pv_dx, ...)	Compute six orthogonal basis fields from composite PV fields.
`compute_strain_basis`(pv_dx_dx, pv_dy_dy)	Compute normal strain deformation basis: ∂²q/∂x² − ∂²q/∂y².
`compute_laplacian_basis`(pv_dx_dx, pv_dy_dy)	Compute Laplacian (diffusion) basis: ∂²q/∂x² + ∂²q/∂y².
`auto_prenorm`(fields[, mask])	Compute per-event pre-normalization constants.
`project_field`(field2d, basis, *[, ...])	Project a tendency field onto the orthogonal basis.
`collect_term_fields`(fields_dict[, include_terms])	Collect and sign-correct PV budget term fields.

RWB Detection

Identifies Rossby Wave Breaking events by detecting overturning (folding) of geopotential-height or PV contours and classifying them as Anticyclonic Wave Breaking (AWB) or Cyclonic Wave Breaking (CWB).

Two classification methods are available, selected via the method argument of detect_rwb_events():

Classification methods
`method`	Algorithm	When to use
`"bay"` (default, recommended)	MATLAB-consistent path-order of max / min latitude intersections across sample meridians.	Recommended for circumpolar-cropped contours and standard event-centred patches. Reproduces Peters & Waugh (1996).
`"tilt"`	Centerline-tilt slope of the overturning envelope. Slope < −0.15 → AWB; slope > +0.15 → CWB; otherwise UNK.	Alternative metric useful for sensitivity analysis. The ±0.15 dead-zone threshold is stored in `TILT_SLOPE_THRESHOLD`.

When analysing RWB on cropped patches the recommended workflow is:

Extract circumpolar contours from the full Northern-Hemisphere Z field using circumpolar_contours(). This finds contours that span ≥ 355° longitude (i.e. truly circumpolar), mimicking the MATLAB blue_circumcontour approach with periodic doubling.
Crop each circumpolar contour to the event-centred patch with crop_contour_to_patch(), returning Δlon/Δlat coordinates.
Pass the full-NH field (field_nh, lat_nh, lon_nh) and patch centre (centre_lat, centre_lon) to detect_rwb_events(), which performs steps 1–2 automatically.

If no full-NH field is supplied, detect_rwb_events() falls back to sampled_longest_contours() on the local patch (legacy mode).

Obtain contours (circumpolar-first or legacy local-patch sampling).
Detect overturning via meridian intersection counting (overturn_x_intervals()).
Classify AWB / CWB using the chosen method (classify_bay() or classify_tilt()).
Apply centroid-x gates to filter spurious detections.
Construct envelope polygons for visualisation.

For 3-D (nlev × nlat × nlon) inputs, reduce_to_2d() extracts a single pressure level or computes an exponential height-weighted vertical average (level_mode="wavg"). The weighting uses weighted_mean_2d() with the atmospheric e-folding height scale.

References: Peters D, Waugh D W (1996) J. Atmos. Sci. 53, 3013–3031. Thorncroft C D, Hoskins B J, McIntyre M E (1993) Q.J.R.M.S. 119, 17–55.

`RWBConfig`([level_mode, try_levels, ...])	Configuration for RWB detection and classification.
`detect_rwb_events`(field2d_patch, x_coords, ...)	Detect and classify all RWB bays in a 2D field.
`circumpolar_contours`(field_nh, lat_nh, lon_nh)	Extract circumpolar contours from a full Northern-Hemisphere field.
`crop_contour_to_patch`(contour, centre_lat, ...)	Crop a circumpolar contour to an event-centred patch, returning relative coordinates (Δlon, Δlat).
`classify_bay`(xline, yline, xa, xb[, n_samp, ...])	Classify a bay as AWB or CWB using path-order sign.
`classify_tilt`(xline, yline, xa, xb[, ...])	Classify a bay as AWB or CWB using centerline tilt slope.
`centerline_tilt`(xm, y_min, y_max)	Slope of the bay centreline (mid-point of envelope).
`overturn_x_intervals`(xline, yline[, ...])	Find x-intervals where contour is overturning (folded).
`envelope_polygon`(xline, yline, xa, xb[, ...])	Construct envelope polygon for an overturning interval.
`sampled_longest_contours`(field2d, x_coords, ...)	Sample the longest contours on a LOCAL patch (legacy helper).
`reduce_to_2d`(arr, levels, level_mode[, ...])	Get 2D slice from 3D array using level index or weighted mean.
`weighted_mean_2d`(arr3d, z3d_m, H_SCALE)	Exponential height-weighted vertical average.
`nearest_level_index`(levels, hpa)	Find nearest pressure level index.
`TILT_SLOPE_THRESHOLD`	Convert a string or number to a floating-point number, if possible.

Composites

Accumulates per-event NPZ patch fields into running sums (grouped by event stage and RWB variant) and exports / loads the composite state via pickle. CompositeState provides composite_mean_3d() and composite_reduce() methods for rapid access to composite-mean fields.

`CompositeState`(globals_map)	Lightweight wrapper around exported composite globals.
`load_composite_state`(path)	Load CompositeState from a pickle file.

Isentropic Interpolation

Interpolates 3-D isobaric fields onto constant-\(\theta\) (isentropic) surfaces using the algorithm of Ziv & Alpert (1994) as implemented in MetPy v1.7.

Algorithm (Newton–Raphson on Poisson’s equation):

Sort pressure levels in descending order (surface → top).
Compute \(\theta = T\,(P_0/p)^\kappa\) at every grid point.
For each target \(\theta^*\), find bounding pressure levels.
Assume \(T\) varies linearly with \(\ln p\) between bounding levels.
Solve for \(p^*\) via Newton–Raphson:

\[f(\ln p) = \theta^* - (a\,\ln p + b)\,P_0^\kappa\, e^{-\kappa\,\ln p}\]
Interpolate all auxiliary fields to the discovered \(p^*\) surface using log-pressure weighting.

Convenience wrappers handle multi-event NPZ data directly.

References: Ziv A, Alpert P (1994). MetPy isentropic_interpolation source.

`isentropic_interpolation`(theta_levels, ...)	Interpolate fields from isobaric to isentropic coordinates.
`isentropic_interpolation_pressure`(...[, ...])	Compute pressure on isentropic surfaces (MetPy algorithm).
`interp_event_fields_to_theta`(theta_3d, ...)	Interpolate one 3-D isobaric field to multiple θ surfaces (vectorized).
`interp_event_field_to_single_theta`(theta_3d, ...)	Interpolate one 3-D field onto a single θ surface → (ny, nx).

Tendency Pipeline

Orchestrates the full per-event computation (Helmholtz-first, v2.0):

Load ERA5 data for the time window (including specific humidity q).
Subtract hourly climatology → anomalies; compute \(\partial T/\partial t\).
Compute all spatial / temporal derivatives, including Emanuel (1987) LHR \(\dot\theta_\text{LHR}\) and \(Q_\text{LHR}\).
Helmholtz on total wind (full NH, spherical Poisson) → \((u_\text{rot}, u_\text{div}, v_\text{rot}, v_\text{div})\).
Load pre-computed Helmholtz climatology (24 monthly NetCDF files, day-resolved 5-D fields since v2.1) → \((\bar u_\text{rot}, \bar u_\text{div}, \bar v_\text{rot}, \bar v_\text{div})\).
Anomaly Helmholtz = total − clim (no separate anomaly solve).
QG omega (SIP, terms A+B) → \(\omega_\text{dry}\).
QG omega (SIP, terms A+B+\(C_\text{LOG20}\) with \(\partial T/\partial t\)) → \(\omega_\text{qg\_moist} = \omega_{A+B+C} - \omega_{A+B}\).
QG omega (SIP, terms A+B+\(C_\text{em}\) with Emanuel LHR) → \(\omega_\text{em\_moist} = \omega_{A+B+C_\text{em}} - \omega_{A+B}\).
Moist / dry decomposition → \(\omega_\text{moist}\), \(\omega_\text{qg\_moist}\), \(\omega_\text{em\_moist}\), and their divergent winds via independent spherical Poisson inversion.
Extract event-centred patches.
Compute 53 PV cross-terms (20 primary + 16 alt-vertical + 16 div dry/moist horizontal + \(Q_\text{LHR}\)) and vertical weighted averages.
Write per-timestep NPZ files.

TendencyComputer is parameterised by event type (blocking / PRP), eliminating code duplication between scripts.

Helper functions

`TendencyConfig`(event_type, data_dir, ...)	Configuration for PV tendency computation.
`TendencyComputer`(config)	Computes PV tendency terms for weather events.

`load_climatology`(clim_path[, engine])	Load climatology, auto-detecting the file layout.
`month_keys_for_window`(base_ts[, hmin, hmax])	Determine which (year, month) files are needed for the time window.
`open_months_ds`(data_dir, var_list, month_keys)	Open multiple months of ERA5 data as a single dataset.
`with_derivs_for_window`(base_ts, cfg, clim_ds)	Open data for base_ts window, compute all bars/anoms/derivatives, Helmholtz decomposition on full NH.
`get_tracked_center`(track_file, track_id, ...)	Look up Lagrangian centre from tracking data.

RWB Classification (Pass 1)

Reads the dh=0 NPZ snapshots produced by TendencyComputer, classifies each event as AWB (Anticyclonic Wave Breaking), CWB (Cyclonic Wave Breaking), or NEUTRAL at multiple pressure levels, and emits a “variant tracksets” PKL.

The classification uses sampled_longest_contours() on Z-field contours followed by overturn_x_intervals() and classify_bay(), with a tilt-based fallback. A multi-level threshold is applied: the event must be classified consistently across at least classify_threshold pressure levels to qualify.

--levels accepts integer hPa values (e.g. 500 400 300 200) or the string wavg to classify on the pre-computed weighted-average 2-D Z field. Use --threshold 1 with a single level or wavg.

The CLI subcommand is pvtend-pipeline classify.

`ClassifyConfig`(npz_dir, output_path, stages, ...)	Configuration for Pass-1 RWB classification.
`ClassifyResult`(stage_all, stage_awb, ...)	Holds RWB variant classification results.

run_pass1(cfg)

Run Pass-1 RWB classification from NPZ files.

Variant-aware Composite Builder (Pass 2)

Accumulates per-event NPZ fields into running sums grouped by event stage (onset / peak / decay), relative hour offset, and RWB variant.

Ten variants are supported:

original — all events (no RWB filter).
AWB_onset, AWB_peak, AWB_decay
CWB_onset, CWB_peak, CWB_decay
NEUTRAL_onset, NEUTRAL_peak, NEUTRAL_decay

CompositeResult provides mean_3d() and reduce_2d() methods for rapid access to composite-mean fields, including exponential height-weighted vertical averages (level_mode="wavg").

The CLI subcommand is pvtend-pipeline composite.

`CompositeConfig`(npz_dir, stages, exclude_file)	Configuration for Pass-2 composite accumulation.
`CompositeResult`(levels, x_rel, y_rel, ...)	Accumulated composite data, supporting original + RWB variants.

build_composites(cfg[, rwb])

Accumulate NPZ fields into variant-aware composites.

I/O

Utilities for loading ERA5 monthly NetCDF files, per-event NPZ patches, and composite-state pickle files.

`io.load_era5_month`(data_dir, year, month, ...)	Load ERA5 monthly data for specified variables.
`io.open_months_dataset`(data_dir, var_list, ...)	Open multiple months of ERA5 data as a single dataset.
`io.load_npz_patch`(path)	Load an event NPZ patch file.
`io.list_npz_patches`(root_dir[, stage, dh])	List all NPZ patch files under a directory.
`io.load_pkl`(path)	Load data from a pickle file.
`io.save_pkl`(data, path)	Save data to a pickle file.

Plotting

Publication-quality visualisation tools for PV tendency analysis.

Basis & coefficient plots

`plotting.plot_four_basis`(phi_int, phi_dx, ...)	Plot the four orthogonal basis fields in a 2×2 grid.
`plotting.plot_basis_with_contours`(basis, *)	Convenience wrapper using an OrthogonalBasisFields object.
`plotting.plot_coefficient_curves`(dh_values, ...)	Plot coefficient time series in a 2×3 panel.
`plotting.plot_field_2d`(field, x, y, *[, ...])	Plot a 2D field on relative coordinates.
`plotting.plot_wind_overlay`(u, v, x, y, *[, ...])	Plot wind vectors (quiver) with optional scalar background.

Composite explorer

plot_var() is a self-contained single-variable composite viewer that loads NPZ events, computes a bootstrap significance mask (N = 1000, 95 % CI by default), and optionally projects onto the orthogonal basis. When projection=True the basis is built from temporally interpolated fields (bi-linear, \(\alpha = 0.75\) by default, i.e. dh − 15 min) using the _next API of compute_orthogonal_basis().

Two layout modes:

projection=False → 2-panel figure: composite mean + hatched bootstrap significance.
projection=True → 6-panel figure: adds 2 × 2 projection rows (INT / PRP / DEF / Residual) with β, αₓ, αy, γ in the subtitle.

Supported options:

Parameter	Description
`data_root`	Path to composite archive (blocking or PRP)
`stage`	`"onset"` / `"peak"` / `"decay"`
`dh`	Lifecycle hour offset (any integer)
`level`	`"wavg"` or int hPa (e.g. 200, 500)
`var_spec`	`str`, `list[str]` (`'-'`-negation), or `callable(event)`
`projection`	`True` for 6-panel, `False` for 2-panel
`use_sig_mask`	If `True` (default), zero out non-significant grid points before projection; set `False` to project the full composite mean.
`mask`	PV anomaly mask specification. `"< 0"` (default) selects q’ < 0 (blocking composites). `"< -2e-7"` for tight single- event masking; `"> 2e-7"` for cyclones. `False`/`None` to disable. A bool 2-D `ndarray` for custom masks.

Note

All bootstrap and composite-mean routines use np.nanmean / np.nanpercentile internally, so events with partial NaN coverage (e.g. high-latitude edge patches) are handled gracefully without corrupting the composite or flipping projection signs.

Helper functions load_events(), get_field(), and bootstrap_sig() are also public so notebooks can do ad-hoc analysis without re-implementing loaders.

`plotting.plot_var`(var_spec, *, data_root[, ...])	Composite + bootstrap plot for any variable(s).
`plotting.load_events`(data_root, stage, dh, *)	Load all event dicts from `data_root/stage/dh=±<dh>/track_*.npz`.
`plotting.get_field`(event, name[, level])	Extract a single named 2-D field from one event dict.
`plotting.bootstrap_sig`(events, var_spec[, ...])	Composite mean + bootstrap significance mask.

Baroclinic tilt overlay

plot_baroclinic_tilt() shows the westward tilt with height characteristic of baroclinic blocking by overlaying upper-level v′ (bold black contours, significant only) on lower-level v′ (blue-red shading). Panels are generated per lifecycle stage.

plotting.plot_baroclinic_tilt(*, data_root)

Plot v′ at two pressure levels to reveal baroclinic tilt.

Sample Data

Bundled idealized PV evolution for quickstart examples and unit tests.

data.load_idealized_pv()

Load the idealized Gaussian PV evolution sample data.