Clarify and complete YSU momentum countergradient integration

[WeiqunZhang]

Summary

The YSU momentum countergradient path in implicit vertical diffusion appears internally inconsistent, but its intended ERF-specific flux convention is unclear. This should be resolved as a contract/integration issue rather than by changing the RHS sign in isolation.

Current behavior

1. The implicit correction is enabled for the wrong YSU implementation

Source/TimeIntegration/ERF_ImplicitPre.H enables momentum countergradient terms only for the old PBLType::YSU:

const bool l_use_ysu_mom_cg =
    solverChoice.turbChoice[level].enable_ysu_countergradient
    && (solverChoice.turbChoice[level].pbl_type == PBLType::YSU);

However:

• ComputeDiffusivityYSU does not populate HGAMU_v or HGAMV_v.
• The old YSU implementation currently rejects unstable conditions, where the nonlocal countergradient would be relevant.
• ComputeDiffusivityYSUNew computes and stores HGAMU_v and HGAMV_v, but PBLType::YSUNew does not enable their use in ERF_ImplicitPre.H.

Thus the implemented momentum countergradient correction appears inactive for the implementation that actually produces the fields.

2. The stored units do not match the documented units

Source/ERF_IndexDefines.H documents:

HGAMU_v // units m/s/m
HGAMV_v // units m/s/m

But Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp stores:

hgamu = brint * u_klo;
hgamv = brint * v_klo;

These have velocity units (m/s), not velocity-gradient units (m/s/m).

The corresponding WRF formulation divides these quantities by PBL height when applying the nonlocal flux. ERF may intentionally use a different equation or normalization, but the current producer, field documentation, and consumer do not establish a consistent contract.

3. The RHS sign needs an ERF-specific derivation

The implicit momentum solvers add the countergradient flux divergence using += at the bottom boundary and in the interior:

• ERF_ImplicitDiff_N.cpp
• ERF_ImplicitDiff_S.cpp
• ERF_ImplicitDiff_T.cpp

The scalar countergradient implementation uses the opposite sign. WRF also applies its momentum countergradient with a leading negative sign, but ERF evolves conservative face-centered momentum and has different density, staggering, and metric factors. Therefore, WRF alone is not sufficient evidence for changing ERF's sign.

The correct sign should be derived from ERF's intended definition of the nonlocal momentum stress. For example, if ERF intends

tau_i3 = -rho K_m (du_i/dz - gamma_i),

then the countergradient contribution to the momentum RHS is -div(rho K_m gamma_i). If HGAMU/HGAMV instead store a signed stress or source term, the current sign may be intentional.

Suggested resolution

1. Document the ERF-specific momentum countergradient equation and sign convention.

2. Decide whether HGAMU_v/HGAMV_v store:

   • raw velocity-scale quantities (m/s), or
   • normalized gradients (1/s).

3. Make the producer, EddyDiff documentation, and all N/S/T consumers use that same definition.

4. Enable the correction for the PBL implementation that actually produces the fields, presumably YSUNew.

5. Re-derive the uniform, stretched, and terrain-fitted RHS terms from the ERF conservative momentum equation, including density and metric factors.

6. Add a manufactured single-column test with prescribed K_m and countergradient values that verifies:

   • the expected direction of momentum transport,
   • column-integrated conservation,
   • agreement across constant-dz, stretched-dz, and terrain-fitted meshes.

Impact

The path appears dormant under the current dispatch, so this issue does not claim a presently reachable reversed-sign runtime failure. However, enabling the existing YSUNew fields without first resolving the sign and normalization contract could produce dimensionally inconsistent or oppositely directed momentum transport.

[hgopalan]

#3486 Will resolve the issue