Notation — symbols & sign conventions

NOTATION — symbols & sign conventions

Linked from every chapter. If the printed equation and the code disagree, the code is wrong. Fix conventions here first.

Coordinates & operators

Rotation & stratification

Dynamical fields

Streamfunction sign. We use $\mathbf u=\hat{\mathbf z}\times\nabla\psi=(-\psi_y,\psi_x)$ and $\zeta=\nabla^2\psi$. In vorticity-form notebooks solve $\nabla^2\psi=\zeta$ (so psi_hat = grid.invert_laplacian(zeta_hat)). Keep this consistent across all notebooks.

Geostrophic, hydrostatic & thermal-wind balance

Geostrophic balance is $\psi=\phi/f$. $u_g=-\phi_y/f,\ v_g=\phi_x/f$ has exactly the streamfunction form above with $\psi=\phi/f$ — geostrophic flow always has a streamfunction, and $\phi/f$ is it. gfdlib.balance works in a single $(y,z)$ cross-section with a closed-form frontal buoyancy field (no discretization): buoyancy_field and its exact $y$-derivative dbdy_frontal, and thermal_wind_u, the exact vertical integral of the thermal-wind relation. The buoyancy front’s sign reverses at $z=H_{trop}$ (an idealized tropopause), so $u_g(y,z)$ has an extremum there at every $y$ — the geostrophic jet core sits where the meridional buoyancy gradient changes sign.

PV inversion & staircases (barotropic, $f$- or $\beta$-plane)

The invertibility principle. Given $q(x,y)$ everywhere, a balance condition ($\zeta=\nabla^2\psi$, i.e. nondivergent flow), and boundary conditions (here: periodic), $\psi$ — and hence $\mathbf u$ — is uniquely determined: psi_hat = grid.invert_laplacian(zeta_hat) with $\zeta=q-\beta y$. gfdlib.pv.gaussian_blob places a periodic-safe PV anomaly at a point; gfdlib.pv.staircase_pv builds an idealized alternating-band PV staircase (a smoothed square wave in $y$). Inverting a staircase shows the jet cores ($u$ extrema) sit exactly at the sharp PV risers between bands, not within the well-mixed plateaus — verified numerically (local $|q_y|$ at each jet core is $\sim$20$\times$ the domain-mean $|q_y|$).

Barotropic instability (Rayleigh-Kuo)

Rayleigh-Kuo equation: $(U-c)(\phi''-k^2\phi)+(\beta-U'')\phi=0$ — derived independently in this repo (not copied from a text) and checked against the classical form at $\beta=0$ (Drazin & Reid). A necessary condition for instability is that $\beta-U''$ changes sign somewhere in the domain (Rayleigh’s inflection-point criterion, generalized to a $\beta$-plane by Kuo). gfdlib.instability.growth_rate solves this as a generalized eigenvalue problem $A\phi=cB\phi$, rewritten $B^{-1}A$ (no SciPy needed) and diagonalized with plain numpy.linalg.eigvals. Verified: a profile with no inflection point gives exactly zero growth at every $k$ tested; the classical $\tanh(y/\delta)$ shear layer’s growth rate peaks at $k\delta\approx0.445$, matching the published value (Michalke 1964, $k\delta\approx0.4446$) to within grid resolution. gfdlib.instability.double_shear_layer builds the periodic nonlinear initial condition; its evolution reuses spectral.Grid and timestep.ifrk4_step unchanged from ch18.

Shallow water (linearized, beta-plane)

Nondimensionalized by $L_R$ (length), $f_0^{-1}$ (time), $H$ (height, so $\hat\eta=\eta/H$, $\hat c=1$), with $\hat f = 1+\hat\beta(y-y_0)$, $\hat\beta=\beta L_R/f_0$:

Linear PV anomaly $q=\zeta-\hat\eta$ (with $\zeta=v_x-u_y$ from the same $\hat u,\hat v$) obeys $q_t=-\hat\beta\hat v$: exactly conserved pointwise on the $f$-plane ($\hat\beta=0$). gfdlib.shallowwater.potential_vorticity computes it; gfdlib.shallowwater.invert_pv inverts $(\nabla^2-1)\eta_{bal}=q$ for the geostrophically balanced height field carrying that PV.

Single-layer QG (barotropic PV + deformation-radius stretching)

QGPV inversion reuses the shallow-water Helmholtz solve. With $L_R$ nondimensionalized to 1, $(\nabla^2-1)\psi=q-\beta y$ is the same operator as the shallow-water PV inversion above — gfdlib.shallowwater. invert_pv computes $\psi$ from the QGPV anomaly directly; no separate QG inversion primitive exists or is needed. The linear Rossby-wave part of the evolution operator changes accordingly: gfdlib.qg.dispersion_omega gives $\omega=-\beta k_x/(k_x^2+k_y^2+1)$, bounded by $|\omega|\le\beta/2$ (attained at $|\mathbf k|=1$) — unlike the unbounded barotropic relation in ch18, the deformation radius screens the lowest wavenumbers. An isolated vortex (gfdlib.pv.gaussian_blob) on this operator sheds a trailing Rossby-wave wake and drifts westward — verified numerically.

Baroclinic instability (Eady; 2-layer/Phillips)

Both growth-rate problems are solved as small eigenvalue problems derived here, not copied from a closed-form formula. gfdlib. baroclinic.eady_growth_rate matches the boundary-value problem’s interior solution to the linearized thermal boundary condition at $z=0,H$; verified against the published benchmark ($\mu_{max}\approx1.61$, $\sigma_{max}\approx0.31$, cutoff $\mu_c\approx2.399$; Vallis/Pedlosky) to within grid resolution. gfdlib.baroclinic.invert_2layer is the 2-layer QGPV Helmholtz solve (closed-form $2\times2$ inverse); rhs_2layer is spectral-space, for use as the nonlinear RHS in gfdlib.timestep.ifrk4_step with hyperviscosity as the (scalar, both-layers-identical) linear operator — the mean-flow and $\beta$ terms stay in the nonlinear part since they couple the two layers ($2\times2$ per wavenumber, not a scalar). Cross-checked against an independently-derived linear eigenvalue: a small-amplitude single-wavenumber perturbation’s measured growth rate converges to the predicted value to $<0.1\%$ once the mode purifies (a coincidental choice of parameters that placed an integer wavenumber exactly at the marginal-stability cutoff was caught and avoided during this check).

A fixed mean shear with no drag is an unlimited energy source. Nonlinearly, rhs_2layer without Ekman drag ($r=0$) was found to grow without bound regardless of how much hyperviscosity was added — energy conservation was verified exact for the pure nonlinear (Jacobian) terms in isolation, so the runaway traced to the mean-flow/$\beta$ terms having no saturating mechanism, not a sign error. Adding bottom drag on the lower layer’s relative vorticity (r>0, the standard ingredient in Phillips’ original 1954 model) produces the expected life cycle — exponential growth, a peak, decay, and equilibration — confirmed by a multi-hundred-time-unit run at the notebook’s actual resolution.

Barotropic/baroclinic mode decomposition (ch. 10). $\psi_{bt}=(\psi_1+\psi_2)/2$, $\psi_{bc}=(\psi_1-\psi_2)/2$ decouple the two-layer PV equations EXACTLY: summing/differencing $q_1=\nabla^2\psi_1+F(\psi_2-\psi_1)$, $q_2=\nabla^2\psi_2+F(\psi_1-\psi_2)$ cancels the $F$ terms in the sum and doubles them in the difference, giving $q_{bt}=\nabla^2\psi_{bt}$ (exactly ch. 7/18’s barotropic PV, no $F$ at all) and $q_{bc}=\nabla^2\psi_{bc}-2F\psi_{bc}$ (exactly ch. 8’s single-layer QGPV, with deformation parameter $2F$). gfdlib.baroclinic. to_modes/from_modes implement the (trivial, linear) change of variables; pv_2layer is the forward operator (exact inverse of invert_2layer). Verified directly, not just asserted from the algebra: to_modes applied to pv_2layer’s output matches the two decoupled single-layer operators to $10^{-8}$. The decoupling holds only for the LINEAR operator — the Jacobian nonlinearity in rhs_2layer still couples the two modes, which is exactly what ch. 10’s notebook demonstrates (a disturbance seeded purely in $q_{bc}$ measurably leaks barotropic KE through nonlinear self-advection alone).

Symmetric, inertial & Kelvin-Helmholtz instability (survey)

Symmetric & inertial instability: one Ertel-PV sign criterion. For a zonal thermal-wind front, gfdlib.symmetric derives (not looks up) the Ertel PV $q=(f-\partial U/\partial y)N^2-f(\partial U/\partial z)^2$ and shows the front is unstable to slantwise (symmetric) displacements iff $qf<0$, i.e. $(1+\mathrm{Ro})<1/\mathrm{Ri}$. Two independently-checkable limits confirm the sign convention: at $\mathrm{Ro}=0$ the marginal curve reduces to the textbook value $\mathrm{Ri}=1$; for $(1+\mathrm{Ro})<0$ (absolute vorticity changes sign) the front is unconditionally unstable at every $\mathrm{Ri}>0$ — pure inertial instability, no stratification dependence. gfdlib.symmetric.classify labels any $(\mathrm{Ro}, > \mathrm{Ri})$ point stable / symmetric / inertial / gravitational ($\mathrm{Ri}<0$) for the Ri-Ro stability map.

Kelvin-Helmholtz: the Taylor-Goldstein equation, and a real numerical pitfall caught along the way. gfdlib.instability. taylor_goldstein_growth_rate solves $(U-c)^2(\phi''-k^2\phi)-U''(U-c) > \phi+N^2\phi=0$ as a natively linear $2n\times2n$ generalized eigenvalue problem (keeping the buoyancy perturbation as an independent unknown rather than eliminating it via division by $(U-c)$, which gives a mathematically equivalent but quadratic-in-$c$ pencil). At $N^2=0$ this reduces exactly (to $10^{-13}$) to gfdlib.instability.growth_rate — a strong internal cross-check, not an external one. For the classic Hazel (1972) profile, growth was initially found NOT to vanish cleanly above the rigorous Miles-Howard threshold $\mathrm{Ri}\geq1/4$ — a hard theorem violation is impossible, so this was root-caused, not dismissed: the spurious growth scaled with grid spacing $\Delta z$ alone (matched $\Delta z$ across different $n$ and domain sizes gave matched spurious magnitude) and shrank as $\Delta z\to0$, confirming genuine 2nd-order truncation error — worst exactly at the delicate $\mathrm{Ri}=1/4$ point, where the equation’s critical-layer indicial roots coalesce, a classically hard point for finite differences (professional codes use spectral methods here). Fixed by tightening the domain/resolution and applying a documented noise floor calibrated to this investigation (_TG_NOISE_FLOOR); the resulting curve decreases smoothly and reaches zero by $J\approx0.4$, correctly reproducing the theorem’s qualitative Ri$\geq1/4$ cutoff (not a knife-edge at exactly $1/4$, an honest limitation of a simple solver on this classically delicate problem). gfdlib.instability.hazel_profile builds $U,N^2$ for this test case.

Rossby-wave ray tracing (sphere, WKB)

Great circles are a theorem here, not a plotting choice. For a background flow $U=\Omega_s a\cos\phi$, both $U$ and the meridional gradient of absolute vorticity scale by $(\Omega+\Omega_s)$ — the same factor — so the stationary wavenumber $K_s^2=2(\Omega+\Omega_s)/(\Omega_s > a^2)$ is independent of latitude, and Hamilton’s ray equations (gfdlib.rossby.ray_rhs, via central differences of dispersion_omega — deliberately not hand-derived symbolically) integrate to exact great circles: great_circle_deviation checks this directly (fit the plane through the ray’s first two sampled points on the unit sphere, then check every later point lies in it) and finds deviations at the $10^{-9}$ level, i.e. floating-point roundoff, not approximation error.

Rayleigh-Benard convection (vertical $x$–$z$ channel)

A genuine bug to avoid re-introducing. With $u=-\psi_z,\ w=\psi_x$, vorticity is $\zeta=w_x-u_z=\nabla^2\psi$ — the same sign as the horizontal streamfunction convention above, not $\nabla^2\psi=-\zeta$. Getting this backward flips the sign of the buoyancy-vorticity feedback and the instability never grows at any $Ra$: verified numerically (see ch14 PR) that the wrong sign gives pure decay for all $Ra$, while the correct sign reproduces subcritical decay ($Ra<Ra_c$) and supercritical growth into finite-amplitude convection ($Ra>Ra_c$) as expected. gfdlib.convection.ChannelGrid (periodic $x$, Dirichlet $z$) and rhs_boussinesq implement the corrected convention; lorenz_rhs is the $(X,Y,Z)$ truncation, sharing gfdlib.timestep.rk4 with every other chapter.

Buoyancy-driven overturning (abyssal recipes; horizontal convection)

Munk’s (1966) abyssal recipe $w\,\partial T/\partial z=\kappa\, > \partial^2T/\partial z^2$ is solved exactly (gfdlib.overturning. abyssal_profile): a linear 2nd-order ODE with constant coefficients, general solution $T=A+Be^{wz/\kappa}$, with $A,B$ fixed by the two BCs via a $2\times2$ numpy.linalg.solve. Large $w/\kappa$ sweeps the bottom value through most of the water column, compressing the transition to the surface value into a thin layer near $z=H$ — the real ocean’s near-uniform abyssal temperature capped by a thin thermocline; verified numerically, not just plotted, against both this limit and the $w\to0$ pure-diffusion (linear) profile.

The overturning cell reuses ch14’s Boussinesq machinery unchanged. gfdlib.overturning.rhs_overturning has the exact same structure as convection.rhs_boussinesq, but replaces the fixed background- stratification advection term ($+\psi_x$, uniform bottom heating) with an explicit, horizontally-varying heating field $Q(x,z)$ (gfdlib.overturning.surface_heating) — warm at one end of the (periodic) domain, cool at the other, concentrated near the surface. This is Rossby’s (1965) “horizontal convection” mechanism: differential heating along a boundary organizes into a single overturning cell (rising at the heated column, sinking at the cooled one) rather than the many small convective rolls of uniform bottom heating — verified directly, not just asserted: the vertical velocity $w=\psi_x$ comes out positive under the heated column and negative under the cooled one in every run tested.

Internal gravity waves (vertical $x$–$z$ plane)

A different $\psi$ than the horizontal chapters. This $(u,w)=(\psi_z,-\psi_x)$ sign convention is for a vertical plane and is unrelated to the horizontal $(x,y)$ streamfunction above — the two never appear in the same notebook. gfdlib.internalwaves reuses spectral.Grid with its $y$-axis standing in for $z$; dispersion_omega/beam_angle give $\omega(k,m)$ and $\theta$; step_leapfrog integrates the $q$ equation and accepts a $z$-dependent $N^2$ directly (no assumption of constant stratification). Nondimensionalized by a reference buoyancy frequency $N_0$ (time) and the domain width (length): $\hat N=N/N_0,\ \hat f=f/N_0,\ \hat\omega=\omega/N_0$.

Eddy transport & mixing (passive tracer in 2D turbulence)

The mean-gradient trick. A passive tracer with an unbounded mean gradient is handled the same way ch16/ch18 handle a mean shear or $\beta y$: split into an imposed linear part and a periodic perturbation $c'$, so gfdlib.mixing.rhs_coupled can evolve $c'$ on the same doubly- periodic spectral.Grid as the vorticity that stirs it, with the mean gradient appearing as an ordinary source term $-\Gamma v$ (the tracer is passive: it’s advected by, but does not feed back on, the vorticity).

$K_{eff}$ is intrinsically $\Gamma$-independent — verified, not assumed. The $c'$ equation is linear in $c'$, so $c'\propto\Gamma$ for a fixed flow realization and $K_{eff}=-\overline{v'c'}/\Gamma$ must be exactly independent of $\Gamma$. Checked numerically: two runs on the identical vorticity history with $\Gamma$ differing by a factor of 2 give $K_{eff}$ agreeing to machine precision (relative difference $10^{-16}$), and a control run with the velocity field identically zero gives $K_{eff}=0$ exactly (no possible eddy flux without a flow). For freely- decaying 2D turbulence, $K_{eff}$ comes out positive (down-gradient) and of the same order as a mixing-length estimate $u_{rms}\times$(domain scale) — a real emergent result of the simulation, not tuned to match.

Wind-driven circulation (Stommel; Munk)

A genuinely new domain shape. Every earlier chapter’s domain was doubly-periodic or a channel with two boundaries; this is the first bounded rectangular basin, four walls. gfdlib.circulation solves the Stommel and Munk gyre problems two independent ways: an exact analytic solution (separating $\psi=f(x)\sin(\pi y)$ reduces the PDE to a linear ODE for $f(x)$ with constant coefficients, solved via its characteristic polynomial’s roots — quadratic formula for Stommel, numpy.roots for Munk’s quartic — and a small numpy.linalg.solve for the boundary constants), and a general 2D finite-difference direct solve (Kronecker- sum operators, numpy.linalg.solve, no SciPy) that works for arbitrary forcing, cross-validated against the analytic solution with clean 2nd-order convergence (gfdlib tests).

Munk’s biharmonic operator is handled by keeping $\zeta=\nabla^2\psi$ as an independent field (the same pattern as ch16/17/19’s coupled states) rather than a hand-derived 4th-order stencil, closed with Thom’s (1933) wall-vorticity formula $\zeta_{wall}=2\psi_{adjacent}/d^2$ (derived by Taylor-expanding $\psi$ from a $\psi=\psi_n=0$ wall) for no-slip on the east/west walls, with free-slip ($\zeta=0$, no correction needed) on the north/south walls — the standard textbook simplification that keeps the problem separable.

A real bug, caught by a convergence test, not an eyeball check. An earlier version enforced no-slip on all four walls (matching physical intuition, but not the separable-solution assumption). It disagreed with the exact analytic solution by a stubborn $\sim30\%$ — but critically, that error did not shrink with grid refinement at any boundary-layer width tested (0.08 to 0.4), immediately distinguishing a genuine discretization error (which must vanish as $\Delta x\to0$) from a boundary-condition mismatch (which does not). The fix — free-slip on north/south, matching the classical Pedlosky/Vallis treatment — restored clean 2nd-order convergence.

Kinematics: strain, vorticity & the Okubo-Weiss parameter

A steady strain + smooth vortex, differentiated by finite differences, not by hand. gfdlib.kinematics.velocity_field is a Lamb-Oseen-like vortex superposed on uniform strain (both individually divergence-free); okubo_weiss gets $S_n,S_s,\zeta$ from central differences of that velocity field, the same safety pattern gfdlib.rossby uses for its ray equations (correctness then rests only on velocity_field itself, not on a hand-differentiated closed form for the smoothed vortex profile). Verified against two independently-derivable limits: with $\Gamma=0$ the flow is exactly linear strain, so central differences are exact (W=4\alpha^2 to floating-point precision, any $h$); at the vortex center with $\alpha=0$ the flow is locally solid-body rotation with zero strain, so $W=-\zeta_{center}^2$ with $\zeta_{center}=\Gamma/(2\pi\sigma^2)$ from the vortex’s small-$r$ Taylor expansion.

Inertial oscillations & the Taylor-Proudman limit

Taylor-Proudman is not a separate solver — it is $b\equiv0$ in Ch. 5’s thermal wind. gfdlib.rotation handles only the inertial-oscillation half (inertial_rhs for RK4 integration, inertial_trajectory for the closed-form circle $w(t)=w_0e^{-ift}$, radius $|w_0|/f$, period $2\pi/f$); the Taylor-Proudman half of ch03 reuses gfdlib.balance.thermal_wind_u directly with the buoyancy contrast swept to zero, since $f\,u_{g,z}=-b_y\equiv0\Rightarrow u_{g,z}=0$ requires no new code, only a new reading of code that already existed. Verified: the numerically RK4-integrated $(x,y,u,v)$ trajectory matches inertial_trajectory to $10^{-6}$; the closed-form trajectory’s radius and speed are exactly constant (a genuine circle) to $10^{-10}$; the thermal-wind shear range is exactly linear in $\Delta b$ (inherited directly from thermal_wind_u’s own linearity — no fitting).

Vertical normal modes (Sturm-Liouville, rigid lid)

The equation is $d/dz[(1/N^2)\,d\Phi/dz]+(1/c^2)\Phi=0$ — no $f_0$ in the eigenvalue problem itself. Derived from scratch (separate the linear hydrostatic equations in $z$, track $w$’s structure through continuity, tie it back to $\Phi'$ through the buoyancy equation) and checked twice: for constant $N$, $\Phi_n=\cos(n\pi z/H)$, $c_n=NH/(n\pi)$ solves it exactly by direct substitution (no $f_0$ anywhere), and mode $n=0$ ($\Phi\equiv$const) is trivially a solution with $c_0=\infty$ — the CONTINUOUS analog of gfdlib.baroclinic’s barotropic mode, no deformation-radius constraint at all.

A real, instructive bug: equal-weight ghost points broke matrix symmetry. A first implementation of gfdlib.stratification. vertical_modes used naive ghost-point differencing for the Neumann BCs; it gave eigenvalues that LOOKED nearly correct (within 0.3% of the constant-$N$ benchmark) but mode SHAPES that were badly, visibly wrong once plotted (values 100x too large in the interior) — because the boundary row’s coefficient came out exactly double its neighboring interior row’s reference to it, making the matrix asymmetric and numpy.linalg.eigh invalid on it (silently: eigh doesn’t check symmetry, it just uses the lower/upper triangle and returns garbage for the rest). Fixed with the standard finite-VOLUME treatment (half-cell control volumes at the two boundaries) symmetrized via a diagonal congruence transform before calling eigh — verified: mode shapes now match the exact constant-$N$ cosines to $10^{-12}$, not just the eigenvalues to a few tenths of a percent. The lesson generalizes: for any self-adjoint eigenvalue problem, checking eigenVALUES against a benchmark is not sufficient — check eigenVECTORS too.

Wave-mean-flow interaction (Stokes drift; Rossby critical layers)

Stokes drift is not a separate force — it’s evaluating $u'$ at the parcel’s true position instead of its mean one. gfdlib.wavemean. stokes_drift is derived from scratch in three lines (displacement $\xi=(A/\omega)\sin(\omega t-kx_0)$, then $\langle\xi\,\partial_xu' > \rangle$) and verified two ways: against direct RK4 particle advection in wave_velocity (Ch. 1’s $d\mathbf x/dt=\mathbf u$ machinery, now genuinely time-dependent) sampled stroboscopically once per period, and by sign (positive for $k,\omega$ same sign — with phase propagation, matching the classical surface-wave result).

Rossby-wave critical layers extend ch. 9’s ray tracing to a Doppler-shifted medium. gfdlib.wavemean.rossby_shear_dispersion adds a background shear $U(y)$ to the barotropic dispersion relation; ray_rhs_shear is its EXACT symbolic ray-equation derivative (safe to hand-differentiate here, unlike ch. 9’s central-difference approach, since the formula is simple). Verified directly by integration, not just asserted from the algebra: $\omega$ stays conserved along the ray to $10^{-6}$ even as $y$ and $l$ both change substantially; the ray approaches but never crosses the analytically-predicted critical layer $y_c=\omega_0/(k\Lambda)$ (linear-shear case), with $|l|\to\infty$ and the meridional group velocity $dy/dt\to0$ as it does — the mathematical signature of a genuine singularity in linear ray theory, not a numerical artifact.

Point vortices (Kirchhoff Hamiltonian system)

Positions are the conjugate pair — phase space is physical space. gfdlib.vortex implements the mutual Biot-Savart RHS, the Hamiltonian and invariants, the symplectic implicit-midpoint step (step_midpoint, fixed-point iteration — leapfrog is unavailable because Kirchhoff’s $H$ is non-separable) and a same-order non-symplectic control (step_heun). Verified against exact solutions: co-rotating pair period $T=2\pi^2d^2/\Gamma$ to $10^{-8}$ and dipole speed $\Gamma/(2\pi d)$ to $10^{-6}$ (both derived by hand, not fitted). The chapter’s punchline is a measured fact, not an assertion: over a 20,000-step chaotic 4-vortex run at the same $\Delta t$, Heun’s energy error drifts secularly to $\sim10^{-3}$ while midpoint’s stays bounded at $\sim5\times10^{-6}$ (tests enforce a 30x separation). Finer structure verified too: linear invariants $Q,P$ are conserved to machine precision by EVERY Runge-Kutta method, while the quadratic $L$ drifts under RK4 ($\sim10^{-8}$) but is exact under implicit midpoint (a Gauss method) up to the fixed-point iteration tolerance — $3\times10^{-14}$ at 12 sweeps, measured.

Wave activity & non-acceleration (barotropic)

The theorem is three one-line derivations chained (linearized PV conservation $\Rightarrow A_t=-\overline{v'q'}$; the Taylor identity $\overline{v'q'}=-\partial_y\overline{u'v'}$; the zonal-mean momentum equation $\bar u_t=-\partial_y\overline{u'v'}$) — see ch. 23. gfdlib.wavemean.wave_activity computes $(\bar u, A, \bar q_y)$ from a vorticity snapshot; $A$ requires $\bar q_y$ bounded away from zero (design experiments with $\beta>\max|\bar U''|$). Verified against the FULL nonlinear solver, not the linear theory it comes from: a small packet on $U_0\cos y$ shear produces $\Delta\bar u$ and $\Delta A$ each $\sim10^{-5}$ that cancel to $<1\%$ of themselves, pointwise in $y$ (0.5% measured at $a=0.05$, $n=128$; the residual is the expected $O(a^3)$ correction, and switching on hyperviscosity breaks the cancellation for the stated physical reason — dissipation destroys pseudomomentum without repaying the mean flow).

Slow manifold & balanced initialization (shallow water)

RMS divergence meters the distance from the slow manifold. gfdlib.shallowwater.divergence is the one-line spectral diagnostic; verified against an analytic case to $10^{-10}$ and exactly zero (to $10^{-12}$) on a geostrophically balanced state. Ch. 24’s measured facts: a fully balanced initialization holds RMS $\delta$ at machine zero for the whole run (the linear slow subspace is exactly invariant on the $f$-plane); partial balance rings with amplitude exactly $\propto(1-\alpha)$ (measured: $\alpha=0.5$ gives precisely half the $\alpha=0$ ringing — linearity, not coincidence); and pointwise PV drift is at machine precision in every run regardless of $\alpha$ — balance is entirely about not exciting the fast modes; the slow (PV) content was never at risk.

Dimensionless numbers

gfdlib.scaling turns the table above into testable functions (rossby_number, ekman_number, froude_number, reynolds_number, burger_number, deformation_radius, coriolis_parameter) plus a coarse three-region classify_regime(Ro, Bu) used by ch04’s regime map: $Ro>1$ unbalanced; $Ro<1,Bu>1$ ($L<L_R$) QG-baroclinic (Ch. 16 territory); $Ro<1,Bu<1$ ($L>L_R$) QG-barotropic (Taylor-Proudman-like, Ch. 3/18 territory). A pedagogical simplification of the classical Ro-Bu regime diagram — real transitions are smooth, not the sharp lines drawn on the map — but the two boundaries ($Ro=1$, $Bu=1$) are exactly the numbers' own definitions, not fitted thresholds.