[Bug]: MRF and YSU counter gradient implementations

[AMLattanzi]

Summary

Countergradient term: metric/normalization/sign errors in the implicit vertical diffusion solvers

Follow-up to the discussion in #3543 and #3544. Auditing the MRF/YSU countergradient
implementation across ERF_ImplicitDiff_{N,S,T}.cpp turns up five distinct issues.
All of them collapse to the correct answer on a uniform, flat grid, which is why the
existing regression tests do not catch them.

Reference commit: 199a6e0 (development, 2026-08-12).

---

1. The governing discretization

In terrain-fitted coordinates with J = h_zeta = dz/dzeta, we have d/dz = (1/h_zeta) d/dzeta, hence

J * dF/dz = dF/dzeta

Since detJ is multiplied through the divergence (as the in-code comments note), the
semi-discrete column equation is

detJ * rho * d(phi)/dt = -(1/dzeta) * ( F_{k+1/2} - F_{k-1/2} )

The MRF/YSU flux is

F = -rho*alpha * ( d(phi)/dz - gamma )   =>   F_diff + F_cg ,   F_cg = +rho*alpha*gamma

met_h_zeta appears in a_tmp/c_tmp only because the diffusive flux contains an inner
vertical derivative (d(phi)/dz ~ dz_inv/h_zeta * delta_phi). The countergradient piece
contains no derivative, so it takes the outer dz_inv (already carried in Fact) and
nothing else — no second dz_inv, no 1/h_zeta.

The sign bookkeeping already in the file confirms the convention:
RHS += Fact * scalar_zflux(klo) is the +F_lo term and RHS -= Fact * F_hi is the -F_hi term.

---

2. Issue list

2.1 _T MRF scalar term carries a spurious 1/met_h_zeta — fixed by #3544

Source/Diffusion/ERF_ImplicitDiff_T.cpp:137, 170

RHS_a(i,j,k) -= Fact * (rhoAlpha_hi * gam_hi / met_h_zeta_hi
                      - rhoAlpha_lo * gam_lo / met_h_zeta_lo);   // wrong

Should be (and this is exactly what #3544 does, matching _N:128):

RHS_a(i,j,k) -= Fact * (rhoAlpha_hi * gam_hi - rhoAlpha_lo * gam_lo);

#3544 is correct and should be merged.

2.2 _T YSU momentum term: 1/met_h_zeta and an extra dz_inv — #3543 as proposed is not correct

Source/Diffusion/ERF_ImplicitDiff_T.cpp:419-422, 455-462

Development currently has the 1/met_h_zeta (issue 2.1 again); #3543 additionally
introduces a dz_inv. Both must go. Dimensional check: Fact * rho*alpha * gamma
with gamma in s^-1 gives kg m^-2 s^-1 = rho*u; an extra dz_inv leaves the term
short by one length.

The correct form is identical to the scalar case:

RHS_a(i,j,k) -= Fact * (rhoAlpha_hi * gam_hi - rhoAlpha_lo * gam_lo);

2.3 HGAMU_v/HGAMV_v are stored un-normalized (root cause of the stray dz_inv)

Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp:751 vs :775-776

hgamt_arr(i,j,0) = HGAMT / pblh;   // K/m   -- gradient-like, correct
...
hgamu_arr(i,j,0) = brint * u_klo;  // m/s   -- NOT divided by pblh
hgamv_arr(i,j,0) = brint * v_klo;

brint is dimensionless, so hgamu has units of m/s. WRF divides by hpbl at the point
of use (phys/module_bl_ysu.F, v4.4.2 L1499):

dsdzu = tem1*(-hgamu(i)/hpbl(i) - ufxpbl(i)*zfacent(i,k)/xkzm(i,k))

This contradicts Source/ERF_IndexDefines.H:219-220, which documents HGAMU_v as
units m/s/m. The dz_inv currently in _N:391,425 and the dz_inv_hi/lo in
_S:407,445 are dimensional stand-ins for 1/h_pbl — wrong by a factor h_pbl/dz,
i.e. O(10-100) too large.

Fix at the source rather than in the solver:

hgamu_arr(i,j,0) = brint * u_klo / pblh;   // s^-1, matches the HGAMT convention
hgamv_arr(i,j,0) = brint * v_klo / pblh;

2.4 YSU momentum countergradient has the wrong sign

_N:391,425, _S:407,445, _T:422,462 all use +=, whereas the scalars use -=.

WRF applies an identical structure to heat and momentum
(f(k) += dtodsd*dsdz, f(k+1) -= dtodsu*dsdz; L1236-1238 for heat, L1499-1502 for
momentum), and ERF stores WRF's raw hgamu including the negative sign of brint.
Momentum should therefore use -=, same as the scalars.

Minor: gfac is 1 whenever stagdir < 2, so multiplying the countergradient term by it
is a no-op inside the guard — suggest dropping it for clarity.

2.5 _S scalar term uses face spacings instead of the cell spacing

Source/Diffusion/ERF_ImplicitDiff_S.cpp:131, 166

In _S, Fact = implicit_fac * dt and the outer divergence factor is the cell spacing
dz_inv, while dz_inv_hi/lo are the face spacings belonging to the inner gradient
(cf. c_tmp = -Fact * rhoAlpha_hi * dz_inv_hi * dz_inv and
RHS += Fact * dz_inv * scalar_zflux). The countergradient term has them swapped:

// current -- wrong on a stretched grid, identical on a uniform one
RHS_a(i,j,k) -= Fact * (rhoAlpha_hi*gam_hi*dz_inv_hi - rhoAlpha_lo*gam_lo*dz_inv_lo);

// correct
RHS_a(i,j,k) -= Fact * dz_inv * (rhoAlpha_hi*gam_hi - rhoAlpha_lo*gam_lo);

2.6 Dispatch gap: the YSU momentum term is unreachable, and reads an unwritten component

Source/TimeIntegration/ERF_ImplicitPre.H:67

const bool l_use_ysu_mom_cg = turbChoice.enable_ysu_countergradient
                           && (turbChoice.pbl_type == PBLType::YSU);

But only ERF_ComputeDiffusivityYSUNew.cpp ever writes HGAMU_v/HGAMV_v;
ERF_ComputeDiffusivityYSU.cpp never touches them. So:

• with pbl_type = YSU, the solver reads a component of eddyDiffs that was never filled;
• with pbl_type = YSUNew, the flag is false and the term is dead.

Either extend the predicate to || pbl_type == PBLType::YSUNew or add an assert — but the
current state should not remain, since it silently reads unwritten data.

---

3. Bottom-face assumption (answering the question raised in #3543)

The lower face at klo carries no countergradient flux. The total surface flux is
supplied by the surface layer model or the Neumann BC, and gamma represents nonlocal
transport interior to the PBL — WRF applies the term only at faces with k < kpbl and
never at the ground. Hence at klo only the upper face contributes:

RHS_a(i,j,klo) -= Fact * rhoAlpha_hi * gam_hi;

This is the same in _N, _S, _T, and for both the scalar and momentum variants, and it
is worth stating in a comment so the asymmetry is not re-flagged later (per @asalmgren's
request in #3543; see also #3486).

Note that HGAM* is stored as a column-constant value that is zero above pbli, so the
face averaging gam_hi = 0.5*(gam_k + gam_{k+1}) produces a half-value ramp at the PBL-top
face. That is fine — arguably better behaved than WRF's hard cutoff — but it is a
deliberate difference worth documenting.

---

4. Suggested verification

Two tests that isolate the metric question from the physics:

1. Metric invariance. Run flat terrain with a uniform vertical offset (z_nd = z + const,
   so h_zeta = 1) against a linearly stretched map (h_zeta = c != 1) at the same physical
   dz. The countergradient tendency must be identical; with the stray 1/h_zeta it scales
   as 1/c.
2. Column conservation. Sum sum_k delta(detJ*rho*phi) from the countergradient term
   alone over a column. The telescoping flux difference must leave only boundary faces, both
   of which are zero for gamma. Any stray metric factor breaks the telescoping and leaks
   theta.

A stretched-grid variant of the second test would also cover 2.5.

---

5. Separate observation (not part of this issue)

ERF's YSU appears to implement entrainment only as an enhanced diffusivity at the PBL top
(ERF_ComputeDiffusivityYSUNew.cpp:1024-1080), not as the explicit
-(w'c')_h (z/h)^3 flux term (hfxpbl/ufxpbl * zfacent in WRF). That explicit entrainment
flux is the term that distinguishes YSU from MRF (Hong, Noh & Dudhia 2006), so it would be
worth confirming the omission is intentional. Filing separately if so.

---

Checklist

• Merge Fix MRF countergradient scaling over terrain #3544 (_T scalar: drop met_h_zeta)
• Rework Fix terrain YSU countergradient scaling #3543: drop both met_h_zeta and the added dz_inv in _T momentum
• Normalize hgamu/hgamv by pblh in ERF_ComputeDiffusivityYSUNew.cpp
• Remove the compensating dz_inv from _N:391,425 and dz_inv_hi/lo from _S:407,445
• Flip the momentum countergradient sign += -> -= in _N, _S, _T
• Fix _S scalar spacing (dz_inv_hi/lo -> dz_inv) at :131,166
• Resolve the YSU / YSUNew dispatch gap in ERF_ImplicitPre.H:67
• Add the metric-invariance and column-conservation regression tests
• Comment the zero-bottom-flux assumption in all six call sites

Steps to reproduce

N/A

Expected behavior

N/A

Git commit hash

199a6e0

Environment

N/A

Relevant logs or output

[hgopalan]

#3486 addressed most of the issues. Do let me know if the PR misses something and I can make the correction.

[AMLattanzi]

@hgopalan #3486 needs to be expanded to correct the implicit diffusion files for the YSU counter gradient term. #3543 was closed because it was wrong. Please add the corrections for that model into #3486 and request a review. #3544 will be merged shortly since it is indeed correct.

[hgopalan]

@AMLattanzi Have pushed the changes. Do you want me to add the regression tests too or was it already provided?

[AMLattanzi]

@hgopalan a metric error still remains. Please correct and add a terrain or stretched grid regression test.

[hgopalan]

Done with adding stretched grid regression test.

[asalmgren]

closed by #3869