SHOC bulk-Richardson PBL height skips interpolation at the first level

[WeiqunZhang]

Severity: low · Category: correctness · Fix order: 135 of 150 — fix this 135th.

_{Filenames are numbered in reverse fix order: 001 = fix last, 150 = fix first. This file is 016.}

Location: Source/PBL/Shoc/ERF_ShocStructure.cpp:156

Based on commit 4d7cf58636893f9 (line numbers refer to that tree).

The defect

What is wrong — The bulk-Richardson PBL-height search special-cases k==1 and assigns pblh = zt(1) instead of interpolating between zt(0) and zt(1); prev_rino = 0 is the correct (and well-defined) Richardson number at k=0, so the interpolation branch is valid there and E3SM/CAM pblintd does interpolate at the first level.

Why it matters

Strongly stable surface layer with dz = 100 m (zt0 = 50 m, zt1 = 150 m) where Ri already reaches 0.9 at k=1: the interpolated PBL height is 50 + (0.3/0.9)*100 = 83 m, but the code reports 150 m. The m_pblh_cc plotfile diagnostic (and the zt < pblh + 200 m gate in ShocTKE) is nearly 2x too deep.

Suggested fix

The search special-cases k==1 and assigns pblh = zt(1) instead of interpolating between zt(0) and zt(1); prev_rino = 0 is the correct and well-defined Richardson number at k=0, so the interpolation branch is valid there — E3SM/CAM pblintd does interpolate at the first level. Remove the special case and let the interpolation handle k==1.

--- a/Source/PBL/Shoc/ERF_ShocStructure.cpp
+++ b/Source/PBL/Shoc/ERF_ShocStructure.cpp
@@ ~line 156: k==1 is special-cased to pblh = zt(1) instead of interpolating between zt(0) and
@@ zt(1). prev_rino = 0 is the correct and well-defined Richardson number at k=0, so the
@@ interpolation branch is valid there — E3SM/CAM pblintd does interpolate at the first level.
-                if (k == 1 || amrex::Math::abs(rino - prev_rino) <= 1.0e-12_rt) {
+                if (amrex::Math::abs(rino - prev_rino) <= 1.0e-12_rt) {

# Removing `k == 1` lets the existing interpolation handle the first level; the degenerate-slope
# guard that remains still protects the division. Verify prev_rino really is initialized to 0
# before the k loop (the claim rests on that).

_{Diff is against 4d7cf58636893f9, hand-written from the current source and not compiled — treat it as the precise intent, not a tested patch.}

Verification evidence

F236 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 140 Real prev_rino = 0.0_rt;, line 156 if (k == 1 || amrex::Math::abs(rino - prev_rino) <= 1.0e-12_rt) { pblh_loc = ztk_agl; } (same at 194). Ri(k=0) is identically 0 (thvk-thv0=0, ztk-zt0=0), so the else-branch interpolation is well-defined at k=1; the short-circuit rounds pblh up to zt(1), overshooting by up to one layer. Unchanged since b7484a7. Reaches the m_pblh_cc plotfile (ShocDriver.cpp:676); the ShocTKE gate (line 236) is also gated on tabs<182 K, so effectively dead there. 700*ustar floor does not mask it for ustar<~0.12.

---

_{Based on commit 4d7cf58636893f9. From an automated audit of Source/: findings were verified against tree 7a597da5 and re-checked at 4d7cf58636893f9 wherever the intervening update touched the code. Audit finding id: F236. Reviewer unit(s): PBL/Shoc-2. Nothing here was compiled or run — the failure scenarios are code reasoning, so the reaching configuration above is the cheapest way to confirm or refute it.}