NEDAS.models.vort3d.core module

Zhu, Smith & Ulrich (2001) minimal 3D tropical cyclone model – dynamical core + moisture/surface-flux/radiation physics + a convective closure.

Sigma-coordinate primitive equations on an f/beta-plane, nz free- atmosphere layers (top-to-bottom) plus one boundary layer fixed at the bottom (nlayers = nz+1 total), interfaces carrying sigma-dot. Finite differences (periodic in x, rigid wall in y), 3rd-order Adams-Bashforth time stepping.

Vertical structure generalized from the paper’s own fixed 3-layer design (nz=2: upper troposphere + lower/mid troposphere, plus the boundary layer) to an arbitrary nz, following the pattern in NEDAS’s qg model (models/qg/python/strat.py): the layer/interface layout is built from arrays (sigma_mid, sigma_int, d_sigma), and the hydrostatic cascade and sigma-dot vertical advection are written as loops over those arrays rather than named per-layer constants. For nz=2, the generated sigma values are forced to match the paper’s own Fig. 1/Table A1 numbers exactly, so this is a strict generalization, not an approximation, of the original 3-layer implementation (re-derived and checked to reduce exactly to the original closed-form 3-layer expressions – see dev log).

Two convective closures are available (convection_scheme):

‘ooyama’ – the paper’s own modified Ooyama (1969) closure (eqs.

16-37). Conceptually defined around exactly “boundary layer + 2 free-atmosphere layers” (cloud-base flux from the boundary layer, entrainment from “the middle layer”, detrainment into “the upper layer”) – not a generic multi-level scheme. Only valid for nz=2.

‘betts_miller’ – a simplified Betts (1986)/Betts-Miller-style column

relaxation, formulated purely in terms of a per-layer reference profile (no named layers), so it works for any nz. Precedented for exactly this purpose: Baik, DeMaria & Raman (1990a) ran a 15-level axisymmetric TC model with the Betts scheme (their Table 2), and ZSU2001 repeatedly compares its own results against theirs (though ZSU2001’s own abstract says “11 levels” for that model – a citation discrepancy against the primary source, not something to match here). See betts_miller_adjustment()’s docstring for the specific (simplified) formulation used here, and the dev log for why this is not literally the same model as ZSU2001 once nz != 2, nor the same model as Baik et al. (their model is axisymmetric radius-height, not this module’s Cartesian x-y channel).

A quick (non-faithful) comparison at nz=14 (15 total layers) against Baik et al.’s reported control-simulation numbers landed in a broadly similar wind-speed range (fluctuating ~50-72 m/s vs. their reported 58 m/s at maturity) but shallower minimum pressure (~955-980 hPa vs. their 923 hPa) and a qualitatively different life cycle (rapid intensification by ~25-50h here vs. their slow 0-48h/rapid 96-144h/steady 144-192h three-stage progression) – expected given the very different vortex/sounding/domain/ geometry, not a discrepancy to chase. See dev log for the full comparison and radius-height cross sections.

See techNotes/models/vort3d.md dev log for the paper’s equations/parameters and running notes on simplifications made here.

NEDAS.models.vort3d.core.theta_from_T_p(T, p)

Potential temperature from temperature T (K) and pressure p (Pa).

NEDAS.models.vort3d.core.sounding_T(sigma)

Far-field environmental temperature (K) at a given sigma, linearly interpolated from the paper’s Appendix A Table A1 (all 5 tabulated sigma positions, not just the nz=2 layer midpoints).

NEDAS.models.vort3d.core.sounding_q(sigma)

Far-field environmental specific humidity (kg/kg) at a given sigma, linearly interpolated from the paper’s Appendix A Table A1.

NEDAS.models.vort3d.core.make_sigma_levels(nz, sigma_boundary_top=0.8888888888888888)

Vertical layout: nz free-atmosphere layers (top-to-bottom) + 1 boundary layer occupying sigma in [sigma_boundary_top, 1] (paper’s own boundary- layer depth, 1/9 in sigma, kept fixed regardless of nz – the paper gives no rule for varying it, and it’s a physically distinct layer, not something that should get thinner/thicker just because nz changes).

For nz=2, returns the paper’s own exact Fig. 1 values (NOT re-derived from a general rule – the paper’s actual interfaces aren’t evenly spaced in sigma, and there’s no way to recover that specific asymmetric choice from a generic formula). For other nz, the free troposphere [0, sigma_boundary_top] is divided into nz EQUAL-sigma-thickness layers – a reasonable default the paper doesn’t specify, since it only ever used nz=2.

Returns:

(nz+1,) layer midpoint sigmas, index 0..nz-1 = free

atmosphere top-to-bottom, index nz = boundary layer.

sigma_int: (nz,) interface sigmas: [0..nz-2] between adjacent

free-atmosphere layers, [nz-1] = sigma_boundary_top (the boundary-layer-top interface).

Return type:

sigma_mid

NEDAS.models.vort3d.core.smith_vortex(r, vm=15.0, rm=120000.0)

Smith et al. (1990) tangential wind profile, eq. (13).

NEDAS.models.vort3d.core.make_grid(nx, ny, dx)

Cartesian (xx, yy) coordinate arrays, shape (ny, nx), centered on the domain (origin at the middle grid point) – where the initial vortex is always placed.

NEDAS.models.vort3d.core.gradient_wind_balance_pstar(r, v_tan, f0, rho0)

Simplified initialization: integrate gradient-wind balance dp/dr = rho*(f*v + v^2/r) radially outward from the vortex center to get the pressure perturbation associated with the initial axisymmetric vortex, instead of solving the paper’s full nonlinear balance equation (their eqs. 14-15, an elliptic PDE via Kurihara & Bender 1980). Documented deviation – see vort3d.md dev log.

NEDAS.models.vort3d.core.qsat(T, p)

Saturation specific humidity via Bolton (1980). T in K, p in Pa; returns kg/kg.

NEDAS.models.vort3d.core.surface_drag_coef(Vb)

Shapiro (1992) neutral drag coefficient, with a velocity cap (documented deviation, see dev log) since the linear formula has no physical ceiling.

NEDAS.models.vort3d.core.random_pressure_field(nx, ny, dx, power_law, seed=None)

Unit-std random pressure perturbation field with a prescribed spectral slope – paired with a geostrophic-wind derivation (Core.__init__) so the background flow starts in mass/geostrophic balance, unlike generating wind directly (the previous random_flow(), which left pstar with no matching perturbation at all – see vort3d.md dev log, 2026-07-21).

class NEDAS.models.vort3d.core.Core(nx=100, ny=100, dx=20000.0, nz=2, Vbg=0.0, Vslope=-3, bg_seed=None, beta=0.0, moist=True, convection_scheme='ooyama', sigma_boundary_top=0.8888888888888888, Vmax=15.0, Rmw=120000.0, vortex_x0=0.0, vortex_y0=-700000.0, u_bkg=0.0, v_bkg=0.0, f0=np.float64(4.988021770261553e-05), theta_offset=0.0, q_offset=0.0)

Bases: object

Standalone (NEDAS-independent) integrator for the vort3d dynamical core: sigma-coordinate primitive equations on an f/beta-plane, nz free-atmosphere layers + 1 boundary layer, 3rd-order Adams-Bashforth time stepping. See the module docstring for the vertical-structure generalization and the two convective closures.

Parameters:

• nx (int) – horizontal grid dimensions (paper: 200x200).

• ny (int) – horizontal grid dimensions (paper: 200x200).

• dx (float) – horizontal grid spacing, m (paper: 20000).

• nz (int) – number of free-atmosphere layers (paper: 2). Total prognostic layers = nz+1 (always +1 boundary layer at the bottom). nz=2 uses the paper’s own exact sigma levels; other nz use equal-sigma-thickness free-tropospheric layers (see make_sigma_levels).

• Vbg (float) – random background-flow wind speed amplitude, m/s (0 = calm, matching the paper’s own experiments).

• Vslope (float) – background-flow kinetic-energy spectrum power law.

• bg_seed (int or None) – RNG seed for the background flow.

• beta (float) – df/dy, Coriolis beta parameter, /m/s (0 = pure f-plane, matching the paper’s own experiments; >0 enables beta-drift).

• moist (bool) – if False, run the dry dynamical core only (no surface fluxes, radiative cooling, condensation, or convection) – valid for any nz.

• convection_scheme (str) – ‘ooyama’ (the paper’s own closure, only valid for nz=2) or ‘betts_miller’ (valid for any nz).

• sigma_boundary_top (float) – sigma at the top of the boundary layer (paper’s own value: 8/9). Only affects nz != 2 (the nz=2 sigma levels are always the paper’s exact Fig. 1 values, regardless of this argument).

• Vmax (float) – initial vortex peak tangential wind speed, m/s (smith_vortex’s own default: 15.0).

• Rmw (float) – initial vortex radius of maximum wind, m (smith_vortex’s own default: 120e3).

• vortex_x0 (float) – initial vortex center, m, relative to the domain center (0,0) – default places it away from the domain center and toward the southern (negative-y) boundary (vortex_y0=-700e3), giving the vortex room to drift poleward/ zonally over a long integration (observed beta-drift + Vbg advection can otherwise run it into a wall if it starts centered – see vort3d.md dev log, 2026-07-22). Configurable (rather than hardcoded to the domain center) so generate_init_ensemble can perturb it per member, mirroring vort2d’s loc_sprd.

• vortex_y0 (float) – initial vortex center, m, relative to the domain center (0,0) – default places it away from the domain center and toward the southern (negative-y) boundary (vortex_y0=-700e3), giving the vortex room to drift poleward/ zonally over a long integration (observed beta-drift + Vbg advection can otherwise run it into a wall if it starts centered – see vort3d.md dev log, 2026-07-22). Configurable (rather than hardcoded to the domain center) so generate_init_ensemble can perturb it per member, mirroring vort2d’s loc_sprd.

• f0 (float) – reference Coriolis parameter, /s, at beta=0 / y=0 (the domain-center latitude) – default 2*7.292e-5*sin(20 deg), matching the paper’s own fixed 20N assumption; beta then adds the y-dependence on top of this reference value.

• theta_offset (float) – domain-uniform (no horizontal gradient) perturbation added ONLY to the boundary layer’s theta/q (K, kg/kg) – all free-atmosphere layers, including the top level, are left exactly at the reference sounding. Meant as a simple ensemble IC-spread mechanism (domain-averaged boundary-layer thermodynamic uncertainty), analogous to vortex_x0/vortex_y0’s position spread – see Vort3DModel.generate_init_ensemble for where these get drawn per member. No accompanying dynamical adjustment needed since a spatially-uniform perturbation has zero horizontal gradient (same reasoning as u_bkg/v_bkg needing no matching pressure term).

• q_offset (float) – domain-uniform (no horizontal gradient) perturbation added ONLY to the boundary layer’s theta/q (K, kg/kg) – all free-atmosphere layers, including the top level, are left exactly at the reference sounding. Meant as a simple ensemble IC-spread mechanism (domain-averaged boundary-layer thermodynamic uncertainty), analogous to vortex_x0/vortex_y0’s position spread – see Vort3DModel.generate_init_ensemble for where these get drawn per member. No accompanying dynamical adjustment needed since a spatially-uniform perturbation has zero horizontal gradient (same reasoning as u_bkg/v_bkg needing no matching pressure term).

• u_bkg (float) – uniform (spatially-constant) steering flow, m/s, added on top of Vbg’s turbulent background flow – the standard simple “steering flow” scheme from the beta-drift/ beta-and-advection TC-motion literature (a constant vector wind, distinct from Vbg’s random field). Unlike Vbg, a uniform flow has zero gradients, so it’s invariant under the model’s (gradient-based) diffusion and doesn’t get sheared apart by the vortex – it actually persists and steers, which a domain-scale random field on a vortex-dominated small domain was observed not to (see vort3d.md dev log, 2026-07-22). No matching pressure perturbation is added: exactly balanced on an f-plane (both sides of geostrophic balance are zero for a spatially uniform field); for beta!=0 there’s a small residual imbalance from f varying with y, neglected here (documented simplification, negligible next to the model’s other approximations).

• v_bkg (float) – uniform (spatially-constant) steering flow, m/s, added on top of Vbg’s turbulent background flow – the standard simple “steering flow” scheme from the beta-drift/ beta-and-advection TC-motion literature (a constant vector wind, distinct from Vbg’s random field). Unlike Vbg, a uniform flow has zero gradients, so it’s invariant under the model’s (gradient-based) diffusion and doesn’t get sheared apart by the vortex – it actually persists and steers, which a domain-scale random field on a vortex-dominated small domain was observed not to (see vort3d.md dev log, 2026-07-22). No matching pressure perturbation is added: exactly balanced on an f-plane (both sides of geostrophic balance are zero for a spatially uniform field); for beta!=0 there’s a small residual imbalance from f varying with y, neglected here (documented simplification, negligible next to the model’s other approximations).

Attributes set after construction: u, v, theta, q (each shape (nz+1, ny, nx)), pstar (shape (ny, nx), column mass p*=ps-p_top). Advance the state in time with step(dt).

ddx(f)

ddy_interior(f)

upwind3_dx(f, u)

upwind3_dy(f, v)

diffuse4(f, k1)

hydrostatic(theta, pstar)

Geopotential at each layer midpoint, via the Arakawa & Suarez (1983) layer-mean-Exner-function scheme the paper adopts (Appendix B), generalized to n=nz+1 layers as a loop instead of literal 3-layer algebra (verified to reduce exactly to the original 3-layer formulas for nz=2 – see dev log for the derivation this was checked against).

For each layer k (0=top free-atm layer .. n-1=boundary layer), bounded by half-levels k (above) and k+1 (below):

Phat[m] = naive Exner at half-level m’s own pressure P[k] = layer-mean Exner (mass-weighted over the layer’s own

bounding half-levels)

theta_hat[m] = interface theta at internal half-level m (weighted

average of the two adjacent layers’ theta, using P)

Phi[n-1] = 0 + cp*theta[n-1]*(Phat[n]-P[n-1]) (lowest layer,

special form: no layer below the boundary layer)

Phi[k] = Phi[k+1] + cp*theta_hat[k+1]*(P[k+1]-P[k]) (k<n-1)

rhs()

Right-hand-side tendencies (du, dv, dtheta, dq, dpstar) for the dynamical core (eqs. 2-8: momentum, hydrostatic PGF, sigma-dot vertical advection, surface fluxes/radiative cooling if moist) – does NOT include diffusion or the convective closure, both applied separately in step(). Returns the tendency tuple used by AB3.

convective_venting(dt)

Modified Ooyama (1969) convective closure (eqs. 16-27, 37). Only called when nz==2 (asserted in __init__), so this references layer indices 0,1,2 directly (=upper troposphere, lower/mid troposphere, boundary layer, exactly the paper’s layer1/layer3/layerb) rather than a generalized loop – the closure itself isn’t a generic multi-level scheme (see module docstring).

betts_miller_adjustment(dt)

Simplified Betts (1986)/Betts-Miller-style convective adjustment, generalized to any nz (see module docstring for why this replaces Ooyama’s closure for nz != 2, and Baik/DeMaria/Raman 1990a,b,1991 as the literature precedent for “Betts scheme + many levels + TC”).

This is a SIMPLIFIED realization, not the full literature scheme (which iteratively adjusts a reference profile shape to conserve column enthalpy and guarantee non-negative precipitation): the reference profile at each free-atmosphere layer is defined as the SATURATED state with the same moist static energy as the boundary layer,

h_ref(k) = cp*T_ref(k) + Phi(k) + Lv*qsat(T_ref(k), p(k)) = h_b

solved iteratively for T_ref(k) (same under-relaxation iteration already used for the Ooyama closure’s downdraft wet-bulb solve – the naive fixed-point form is not contractive here either, for the same reason: |Lv/cp * dqsat/dT| typically exceeds 1). This represents “what the column would look like if fully neutralized by deep convection sourced from the boundary layer” – a reasonable, defensible simplification of a true moist-adiabat reference profile, not an exact one. theta and q in each free- atmosphere layer then relax toward theta_ref=T_ref*(p0/p)^kappa and q_ref=qsat(T_ref,p) over tau_BM, wherever the same boundary-layer-instability trigger used by Ooyama’s closure (h_b > max(hstar over free-atm layers)) is satisfied. The boundary layer itself is not relaxed (matches the classic scheme, where the boundary layer is the moisture/energy SOURCE for the reference profile, not itself adjusted).

Does not explicitly conserve column moist static energy (theta and q relax independently, and Phi isn’t updated mid-relaxation) – documented approximation, consistent with the many other documented simplifications already in this model.

step(dt)

Advance the model state by one time step dt (seconds), in place: 3rd-order Adams-Bashforth for the dynamics (rhs(), automatically dropping to lower-order AB for the first two calls while _prev_tendencies fills up), then 4th-order horizontal diffusion as a separate explicit-Euler correction (split from the AB3-integrated tendency for stability, see dev log), then – if moist – condensation and the selected convective closure.

diagnostics()

Per-layer summary statistics (max wind speed, theta range, q range, all as plain lists indexed by layer 0..nz-1=free atmosphere then nz=boundary) plus the domain pstar range – a quick sanity/ progress check, not a substitute for the full state.