pvtend.solve_qg_omega_sip

pvtend.solve_qg_omega_sip(u: ndarray, v: ndarray, t: ndarray, lat: ndarray, lon: ndarray, plevs_pa: ndarray, *, center_lat: float | None = None, omega_b: ndarray | None = None, t_full: ndarray | None = None, rhs_c: ndarray | None = None, phi_3d: ndarray | None = None, bc_top: ndarray | float | None = None, bc_bot: ndarray | float | None = None, bc_lateral: ndarray | None = None, maxit: int = 300, alpha: float = 0.93, resmax: float = 1e-14, w_physical_max: float | None = None) → tuple[ndarray, dict][source]

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

Matches the Li & O’Gorman (2020) solver from dante831/QG-omega/source/SIP_inversion.m.

Solves the σ-in-Laplacian form:

\[\nabla^2(\sigma\omega) + f_0^2 \frac{\partial^2\omega}{\partial p^2} = A + B + C\]

where A is the Dostalek (2017) spherical Q-vector divergence, B = f₀β∂vg/∂p (direct form), and C = -(κ/p)∇²J (diabatic).

The SIP stencil uses local σ(k,j,i) at neighbor points in the horizontal coefficients (AW, AE, AS, AN), and no σ in the vertical coefficients (AB, AT). The meridional coefficients (AS, AN) are discretised in conservative cos(φ_{j±½}) flux-divergence form (Lynch 1989; cf. xinvert), replacing the previous tan φ / Δφ coefficient form which was ill-conditioned at high latitudes. The polar metric vanishes naturally as cos(φ_{j±½}) → 0 at φ → ±90°.

Polar boundary handling (full-NH ring, periodic_lon=True): when the lateral boundary row sits at φ = ±90°, only the m = 0 (zonal-mean) spectral mode is geometrically well defined, so the polar Dirichlet row is replaced by the longitudinal mean of omega_b (or bc_lateral). This matches the spherical-FFT Helmholtz solver (where m ≠ 0 modes vanish at the pole through the cos⁻²φ eigenvalue scaling) and the SH-based _solve_chi_nh divergent-wind inversion (NH→global parity mirror with closed pole).

Boundary conditions: Faces can be individually controlled:

  • omega_b: prescribe ERA5 ω on all Dirichlet faces (default mode).

  • bc_top/bc_bot: override top/bottom faces (e.g. 0.0 for dry solve).

  • bc_lateral: override N/S lateral faces (e.g. ω_bar climatology).

Diabatic heating (term C): When rhs_c is provided (pre-computed via _compute_diabatic_rhs_log20() or _compute_diabatic_rhs_emanuel()), it is added directly to the A+B RHS before solving.

Geostrophic temperature: When phi_3d (geopotential) is provided, thermal-wind angular gradients of T_geo are used for the Q-vector, matching LOG20 GEO_T=true.

Parameters:

  • u – Geostrophic zonal wind [m s⁻¹], (nlev, nlat, nlon).

  • v – Geostrophic meridional wind [m s⁻¹], (nlev, nlat, nlon).

  • t – Temperature [K], (nlev, nlat, nlon).

  • lat – Ascending latitude [degrees], (nlat,).

  • lon – Longitude [degrees], (nlon,).

  • plevs_pa – Pressure [Pa], ascending, (nlev,).

  • center_lat – Latitude for constant f₀ and β₀.

  • omega_b – Boundary omega (nlev, nlat, nlon) (None → 0).

  • t_full – Full physical temperature [K] for σ₀ and σ_3d.

  • rhs_c – Pre-computed diabatic forcing C (nlev, nlat, nlon), -(κ/p)∇²J. Added directly to A+B RHS (σ-in-Laplacian).

  • phi_3d – Geopotential Φ [m² s⁻²], (nlev, nlat, nlon). When provided, T_geo from hydrostatic relation replaces actual T for Q-vector ∇T gradients.

  • bc_top – Top boundary (k=0). 2D array or scalar.

  • bc_bot – Bottom boundary (k=-1). 2D array or scalar.

  • bc_lateral – Lateral boundary omega (nlev, nlat, nlon). Overrides N/S faces from omega_b.

  • maxit – Max SIP iterations (default 300).

  • alpha – SIP relaxation parameter (default 0.93).

  • resmax – Convergence threshold (default 1e-14).

  • w_physical_max – Optional hard clip on abs(omega) at the 300 hPa level in [Pa s⁻¹]. When None (the default), no clipping is applied inside the solver; instead, blowup timestamps are flagged post-hoc by pvtend.blowup.scan_omega_blowup(), which applies the canonical ±5 Pa/s hard cutoff at 300 hPa and emits a CSV of offending timestamps for downstream exclusion. Earlier versions of pvtend silently clipped ω here at ±5 Pa/s, which corrupted the affected events; pass an explicit float only for legacy reproducibility.

Returns:

(omega, info) where omega is (nlev, nlat, nlon) QG vertical velocity [Pa s⁻¹] and info is a dict with 'iters', 'final_residual', 'numba', 'solve_time', 'terms'.

References

Stone H L (1968). Iterative solution of implicit approximations of multidimensional partial differential equations. SIAM J. Numer. Anal., 5, 530–558.

Li L, O’Gorman P A (2020). Response of vertical velocities in extratropical precipitation extremes to climate change. J. Climate.