Author: Maria Chinita and Claude Opus 4.7 Date: 2026-05-23
PR: maint-3.0-scm-mchinita-shoc_calculate_lscale → maint-3.0-scm (PAESCAL-SciDAC5/E3SM-fork)
Commits: 4b1fd1d6ee, 5be2ebd3ce, d8acd77cf5
Three related changes to SHOC, added via a new module shoc_lscale_mod.F90 plus edits to shoc.F90 and shoc_intr.F90. Two are pure diagnostics (bit-for-bit); one is an answer-changing fix for nadv > 1. All validated in SCM on DYCOMS RF01, in two configurations: V1 & V4 use the standard L80 grid, 4 h, 30-min host step (6 macmic substeps); V2 & V3 use the dz=10 m grid (sigma_coord_with_ptop_dz010m_zmax3000m), dt=60 s, 6 h, shoc_outer_loop=true.
/global/cfs/projectdirs/m4359/mchinita/shoc_lscale_thvfix_shocmix_RF01/ — runs/ (simulation outputs), figures/, scripts/, docs/.runs/. PDF originals under CFS figures/.| # | Change | Where | Type |
|---|---|---|---|
| 1 | CLUBB nonlocal length-scale diagnostic shoc_compute_lscale → history LSCALE_SHOC/_UP/_DOWN |
new shoc_lscale_mod.F90, called from shoc_main |
Diagnostic (BFB) |
| 2 | Refresh shoc_thv each nadv substep (host thv goes stale as thetal/qw evolve); built from diagnosed T so it keeps the Exner factor the crude thv_now dropped, matching shoc_intr's thv formula |
shoc.F90 (shoc_main) |
Answer-changing for nadv > 1 |
| 3 | Decompose SHOC_MIX into its 3 component length scales → history SHOC_MIX_SURF/LINF/STRAT |
new compute_shoc_mix_terms in shoc.F90 |
Diagnostic (BFB) |
New h0 fields (add_default): LSCALE_SHOC, LSCALE_UP_SHOC, LSCALE_DOWN_SHOC, SHOC_MIX_SURF, SHOC_MIX_LINF, SHOC_MIX_STRAT.
LSCALE_SHOC closely reproduces the host-side calculate_lscale — the L, L_up, and L_down profiles overlay, with small differences that concentrate at the inversion.
The small residual — a <1% bias through the bulk plus an isolated L_down spike at the inversion — comes from how each routine builds the virtual potential temperature thvm: the host calculate_lscale uses CLUBB's linearized thvm, whereas shoc_compute_lscale uses SHOC's exact-product θv. A controlled experiment that forced calculate_lscale to use the exact-product form collapsed the difference (L total 10.55→0.21 m; L_down 57.68→0.38 m), pinning the cause to the thvm construction — not the TKE path, thlm, or the Exner factor. (Override archived at scripts/lscale_thvmtest_SourceMods/lscale_mod.F90.)
Run (CFS archive): runs/V1_lscale_validation/.
L, L_up, L_down profiles at t=4 h — host reference calculate_lscale (solid) vs SHOC-side shoc_compute_lscale (dashed) lie on top of each other (close agreement).
The experiment that proves the cause. 2×2 difference curtains (LSCALE_SHOC − LSCALE; top L total, bottom L_down). LEFT = original (host CLUBB-linearized thvm): max|Δ| 10.55 / 57.68 m. RIGHT = host forced to SHOC's exact-product thvm form: 0.21 / 0.38 m (featureless). Changing only the thvm construction removes the difference.
Difference curtains (as-shipped): total-L (m and %), and split by up/down component — <1% in the bulk; the only sizable departure is the isolated L_down inversion spike.
shoc_thv == reference thv.Runs (CFS archive): runs/V2V3_shocthv/{mac_only,turb_standalone}/; baselines runs/V2V3_baseline_brigaudet/. 4-curve profiles at t=6 h: baseline mac_only, baseline turb_standalone (stale θv), shoc_thv mac_only, shoc_thv turb_standalone.
u(z) — the clearest: the stale-θv turb curve bows away while the shoc_thv turb curve collapses back onto mac_only.
TKE(z) — same story: the stale curve spreads to ~770 m; the fix recovers the ~635 m structure.
Scalars (T, qt, qv, qc) all four curves overlay — the bug doesn't register in the scalars (well-mixed, no surface fluxes); see the V2V3_shocthv/ portal folder (*.shocthv_4curve.png).
Both fixes refresh thv each substep; they differ only at the cloud-band edges, ~1% of the staleness bug (u ~6e-3 m/s, T ~2e-5 K). Both validate equally; shoc_thv is the physically-correct refinement. Runs (CFS archive): runs/V2V3_thv_now/ vs runs/V2V3_shocthv/.
Δ(shoc_thv − thv_now) for u/T/tke vs height — tiny, confined to the cloud-band edges (where the Exner term acts).
Per-variable 4-curve overlays (mac/turb × thv_now/shoc_thv) all coincide → the two fixes give equivalent results; see the V2V3_shocthv/ portal folder (*.shocthv_vs_thvnow_4curve.png).
The quadrature reconstruction (1/SURF² + 1/LINF² + 1/STRAT²)^(−½) reproduces SHOC_MIX to < 1% (residual = time-averaging nonlinearity); components physically sensible (wall binds near the surface, stratification aloft). Run (CFS archive): runs/V4_shocmix_terms/.
Decomposition + reconstruction (2×2): top row SHOC_MIX with its three component lengths (quadrature; smallest wins) + reconstruction overlay; bottom row LSCALE with its up/down parts (geometric mean) + its reconstruction.
SHOC_MIX decomposition + quadrature reconstruction matching the model's SHOC_MIX to <1%.
pblintd internal thv uses the crude no-Exner reconstruction (same as thv_now). Its thv only sets pblh, which feeds the solution solely through the eddy_diffusivities stable branch gated at tabs(nlev) < 182 K (inactive for DYCOMS/warm cases) → zero solution impact. Left unchanged by decision; documented as a known convention difference.
nadv > 1; validated in SCM, full-model/climate impact not yet characterized.SHOC_MIX itself is unchanged.