Donelan and Variable Rougness MOST Fixes

[AMLattanzi]

Summary

The PR is about the sea-surface roughness pathway. It fixes six correctness
defects, replaces the iteration scheme in two MOST solvers, removes a redundant
flux type, and consolidates five roughness laws into one place so the bounds and
the smooth-flow limit are written once instead of five times.

---

1. Correctness fixes

1.1 Wave-coupled roughness used turbulent, not molecular, viscosity

The aerodynamically-smooth term in the wave-coupled law read

0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar

EddyDiff::Mom_v is assigned as rho * nu_turb
(ERF_ComputeTurbulentViscosity.cpp:148) — dynamic turbulent viscosity,
O(1) kg m⁻¹ s⁻¹. Smith (1988) calls for molecular kinematic viscosity,
ν ≈ 1.5 × 10⁻⁵ m² s⁻¹. The expression was five orders of magnitude too large and
dimensionally wrong (kg m⁻², not m): at u* = 0.3 m/s it evaluated to ≈ 0.37 m,
which pinned z₀ at the z0_max clamp and made the Taylor–Yelland wave-steepness
term irrelevant. The wave-coupled option was, in effect, not doing wave coupling.

Now air_viscosity(tm_arr(i,j,k)) — the Andreas (1989) kinematic ν, 1.46 × 10⁻⁵
m² s⁻¹ at 15 °C. Affected adiabatic_wave_coupled, surface_flux_wave_coupled,
and surface_temp_wave_coupled.

1.2 Charnock smooth-flow term was missing its coefficient

// before
z0 += air_viscosity(tm_arr(i,j,k)) / std::max(ustar, amrex::Real(0.05));
// after
z0 += amrex::Real(0.11) * nu / std::max(ustar, amrex::Real(0.01));

The 0.11 of z₀ = 0.11 ν/u* (Smith 1988; ECMWF IFS; WRF) had been dropped, making
the smooth-flow contribution ~9× too large. The u* floor also moves 0.05 → 0.01
to match the Donelan branch and WRF.

1.3 Moisture flux used a stale roughness

In surface_flux_charnock, surface_flux_mod_charnock, surface_flux_donelan
and surface_flux_wave_coupled, momentum and heat used the freshly updated local
z0 while moisture used the stored z0_arr:

- (std::log(zref / z0_arr(i,j,k)) - psi_h)
+ (std::log(zref / z0) - psi_h)

z0_arr is only written at the end of iterate_flux, so the latent heat flux
lagged the sensible heat flux by a full surface-layer update.

1.4 Modified Charnock was unbounded

z0 = exp[ (2.7 u* − 1.8/b) / (u* + 0.17/b) ],   b = (1/30) ln(1260/d)

has an asymptote of e^2.7 = 14.9 m and previously had no clamp at all. In
shallow water at storm winds the inner z₀ fixed point fails outright — either
iter_z exceeds max_iters and the device-side AMREX_ALWAYS_ASSERT aborts, or
z₀ climbs above the reference height, ln(z_ref/z₀) goes negative, and u*, θ*
and every derived flux are silently corrupted:

z_ref	U₁₀	depth	unbounded result
10 m	20 m/s	5 m	MAXITER (101) → abort
10 m	30 m/s	5 m	z₀ = 12.7 m > z_ref → C < 0
10 m	30 m/s	30 m	MAXITER (101)
10 m	40 m/s	30 m	z₀ = 13.2 m > z_ref
20 m	40 m/s	5 m	64 iters → z₀ = 7.3 m

Now bounded to [1.27 × 10⁻⁷, 1.0 × 10⁻²] m. With that bound every case above
converges in ≤ 24 iterations against max_iters = 100.

1.5 bulk_coeff never produced a roughness length

The bulk-coefficient flux type prescribes C_d/C_h/C_q and computes u*, θ*, q*
from the resulting fluxes, but never wrote z_0. Downstream consumers — the YSU
and YSU-New PBL schemes, 2D plotfiles, checkpoints — saw whatever z0_const was
set at initialization. Compute_roughness now back-solves the log law from the
diagnosed u* and L:

z0 = z_ref / exp( κ U / u*  +  ψ_m(ζ) )

so the reported roughness is consistent with the stress actually applied. It
carries a lower bound only (1 × 10⁻⁷ m) — no ocean cap, because the bulk_coeff
path runs over land as well, and a floor to stop the exponential underflowing to
exactly zero when u* is small. ERF_ComputeDiffusivityYSU.cpp:133 divides by z₀
without a guard, so a zero there would be fatal.

1.6 COARE 3.0 was not COARE 3.0

erf.most.charnock_viscosity now defaults to true (was false). COARE 3.0
(Fairall et al. 2003) is defined with z₀ = α u*²/g + 0.11 ν/u*; the smooth-flow
term is part of the parameterization, not an option, and WRF's module_sf_mynn.F
applies it unconditionally. With the old default, roughness_type_sea = coare3.0
silently ran a variant that isn't COARE.

---

2. Rewritten iteration in surface_temp_donelan and surface_temp_wave_coupled

These two structs iterated directly on u* with Businger–Dyer stability functions
and an ad-hoc bisection guard for Obukhov-length sign flips:

Olen = -u*^3 * thv / (kappa * g * tflux);
if (sign change && |L| + |L_old| < 1.0) Olen = 0.5*(Olen + Oleno);
while (|u*_new - u*_old| > 1e-5)

That heuristic is fragile near neutral, where ⟨w'θ_v'⟩ → 0 drives L → ±∞ and the
iterate chatters. They now use the bulk-Richardson fixed point that
surface_temp, surface_temp_charnock and surface_temp_mod_charnock already
used, so all five specified-temperature solvers share one scheme.

1. First call assumes near-neutral (L = 10³ m) rather than guessing u*.

2. Optional Beljaars (1995) gustiness, U²_eff = |U|² + (β w*)², β = 1.2, with the
   buoyancy flux lagged one update as in WRF.

3. Bulk Richardson number over the reference height, moisture carried through θ_v,
   clipped to [−4, 4].

4. Outer fixed point on ζ using the exact MOST identity

   $$\mathrm{Ri_b} = \zeta,\frac{C-\psi_h}{(C-\psi_m)^2} \quad\Longleftrightarrow\quad \zeta = \mathrm{Ri_b},\frac{(C-\psi_m)^2}{C-\psi_h},\qquad C = \ln(z_{\rm ref}/z_0)$$

   under-relaxed at α = 1/2, with max(C−ψ, 1) guards.

5. New inner fixed point on z₀, since every roughness law depends on u* which
   itself depends on ln(z_ref/z₀). Previously z₀ lagged by one step.

Ri_b is built entirely from known quantities, so the outer iteration has a single
unknown and no singularity at neutral. Operational schemes (WRF sf_sfclayrev,
MYNN) invert Ri_b → ζ the same way, though through regression fits; solving the
relation directly keeps it exact for whatever ψ functions are supplied.

Both structs also switch from Businger–Dyer to the Jiménez et al. (2012) ψ
functions, matching the other surface_temp_* variants. Stable branch follows
Cheng & Brutsaert (2005) with 6.1/2.5 and 5.3/1.1; unstable blends Kansas
(16 ζ) with the Grachev et al. (2000) convective forms (10 ζ, 34 ζ) through
ζ²/(1+ζ²). All coefficients match the published values.

---

3. Removal of the donelan flux type

erf.surface_layer.flux_type = donelan and the donelan_flux functor are gone.
The functor applied a wind-dependent drag coefficient directly:

Cd = 0.001                     U <= 5 m/s
Cd = 0.001 + 7e-5 (U - 5)      5 < U < 25
Cd = 0.0024                    U >= 25
Ch = 0.0012 (constant)
Cq = 0      (no latent heat flux at all)

The ramp is a Large & Pond (1981) style fit and the 0.0024 ceiling encodes the
Donelan et al. (2004) drag saturation, so the shape was defensible. Three
problems made it not worth keeping:

• compute_q_flux returned zero unconditionally — moist runs on this path lost
  all surface latent heat.
• The coefficients were hard-coded; bulk_coeff_flux already offers the same
  functional form through erf.most.Cd/Ch/Cq.
• It bypassed MOST entirely, so drag depended on wind speed alone with no
  stability correction.

Donelan drag is now reached through MOST with
erf.most.roughness_type_sea = donelan, which recovers a wind-dependent C_d
implicitly through z₀(u*) and additionally gets stability corrections and a real
moisture flux. No inputs file in Exec/ or Tests/ used the removed flux type.
The parser aborts on unknown strings, so a stale external deck fails loudly.

---

4. Roughness bounds

Every law now returns through an explicit clamp. What changed:

law	before	after	binds at
Charnock	unbounded	[1.27e-7, 2.85e-3]	U₁₀ ≈ 24.5 m/s
COARE 3.0	[1.27e-7, 2.85e-3]	unchanged	U₁₀ ≈ 24.8 m/s
Donelan	[1.27e-7, 2.85e-3]	unchanged	U₁₀ ≈ 31.4 m/s
modified Charnock	unbounded	[1.27e-7, 1.0e-2]	depth-dependent
wave-coupled	[1.0e-6, 1.0e-1]	[1.27e-7, 1.0e-2]	H_s/L_p ≈ 0.07
log-law inversion	n/a (new)	[1.0e-7, ∞)	floor only

The 1.27 × 10⁻⁷ / 2.85 × 10⁻³ pair is WRF's open-ocean bound
(ZNT = MIN(2.85e-3, MAX(ZNT, 1.27e-7)) in module_sf_sfclayrev.F90), inherited
via Davis et al. (2008). It is best read as a drag ceiling rather than a length:

z₀	C_d10 = (κ/ln(10/z₀))²
1.27 × 10⁻⁷ m	5.09 × 10⁻⁴
2.85 × 10⁻³ m	2.52 × 10⁻³ — the Donelan et al. (2004) saturation value
1.00 × 10⁻² m	3.52 × 10⁻³

Charnock, COARE 3.0 and Donelan keep the WRF value. The two laws whose purpose is
to be rougher than deep water — modified Charnock (shallow water,
Jiménez & Dudhia 2018) and wave-coupled (steep young seas, Taylor & Yelland 2001)
— get 1 × 10⁻² m instead. The deep-water cap would otherwise switch them off in
exactly the regime they exist for: modified Charnock in 5 m of water saturated at
U₁₀ ≈ 10 m/s, and Taylor–Yelland crossed 2.85 × 10⁻³ at a steepness of
H_s/L_p ≈ 0.05 (H_s = 4 m, L_p = 80 m gives 6.7 × 10⁻³ m).

The lower bound is nearly inert with the smooth-flow term active: z₀ ≥ 0.11ν/u*
only falls below 1.27 × 10⁻⁷ at u* > 12 m/s. It matters when visc = false in
near-calm, and it keeps ln(z_ref/z₀) finite.

---

5. Consolidation

ERF_MOSTRoughness.H is deleted and its contents, plus most_data,
similarity_funs and air_viscosity from the top of ERF_MOSTStress.H, move to
the new ERF_MOSTUtils.H. The five roughness laws become named functions —
Charnock_roughness, Mod_Charnock_roughness, COARE3_roughness,
Donelan_roughness, WaveCoupled_roughness — plus Compute_roughness for the
log-law inversion. Previously each law was written inline in three places
(adiabatic / specified-flux / specified-temperature), which is how 1.1, 1.2 and
1.3 came to differ between call sites.

ERF_MOSTStress.H drops from ~2450 to ~2130 lines and now holds only the MOST
functors. Make.package registers the new header; CMake picks it up through the
recursive glob at CMakeLists.txt:724.

---

6. Behavior changes reviewers should expect

These are intended, but they move results:

• charnock_viscosity renamed to a general smooth_flow_viscosity option that changes default false → true. Affects all sea iterators.
• Charnock and modified Charnock are now clamped where they were not.
  Charnock changes above U₁₀ ≈ 25 m/s; modified Charnock above the depth-dependent
  threshold. MOST_test_suite sea1 is modified Charnock.
• Wave-coupled cap tightened 0.1 → 0.01 m. Relevant to
  Exec/CanonicalTests/ABL/WW3_coupling.
• Charnock smooth-flow term is ~9× smaller (the missing 0.11), which shifts
  low-wind z₀ wherever visc is on.
• Wave-coupled z₀ changes by orders of magnitude — the eddy-viscosity bug had
  it pegged at the clamp; it now tracks wave steepness as intended.
• flux_type = donelan is a hard error. External decks using it must switch
  to bulk_coeff, or to MOST with roughness_type_sea = donelan.
• The max(C - ψ, 1) limiter now also sets the final diagnosed u*, t* and q* in the constant-roughness surface_temp solver and in the EB path (Source/EB/ERF_EBMOSTStress.H). Previously the guarded denominators fed only the ζ fixed-point update while the stored coefficients used the raw C - ψ. That left the stored state inconsistent with the converged iteration whenever the guard bound.
• The guarded denominator is now used inside the inner z₀ iteration of all four surface_temp_* variable-roughness solvers, not just for the outer ζ update. Without it the converged z₀ was the roughness-law fixed point of an unguarded u* while the stored u* came from the guarded one, so the stored pair (z₀, u*) did not satisfy the roughness law when the limiter bound.
• erf.most.charnock_viscosity is renamed erf.most.smooth_flow_viscosity and is now read for every sea roughness law, not only charnock. It is queried once before the roughness_type_sea branch and plumbed into all twelve sea functor constructors, so explicit roughness_type_sea = coare3.0, donelan and wave_coupled now honour it as the documentation implies; previously the flag reached only the three *_charnock constructors and the others silently fell through to most_data::visc{true}. The old input name is no longer queried, so external decks setting erf.most.charnock_viscosity are ignored without warning and must be updated. Note that the flag has no effect for modified_charnock: Mod_Charnock_roughness carries its smooth-flow limit intrinsically and takes no visc argument.
• erf.most.modified_charnock_depth is clamped to [10, 100] m. Cnk_b = (1/30) ln(1260/d) vanishes at d = 1260 m and turns negative beyond it, making 1.8/Cnk_b and 0.17/Cnk_b singular and letting the denominator u* + 0.17/Cnk_b cross zero for positive u*; the trailing min/max clamp on z₀ does not sanitize the resulting NaN. The bounds are the fit range of Jiménez & Dudhia (2018), which keeps Cnk_b ∈ [0.084, 0.161]. most_data::Cnk_b also gains an in-class initializer, removing an uninitialized read for functors that never set a depth.

At typical ABL forcing (u* ≈ 0.3–0.5 m/s) none of the new clamps bind, so the
canonical ABL cases are expected to be unaffected apart from the
charnock_viscosity default — but that is worth confirming rather than assuming.

---

7. Known remaining items

Not defects, but worth recording:

• cnk_visc is plumbed only into the three *_charnock constructors. The
  donelan, mod-charnock and wave-coupled structs fall through to
  most_data::visc{true}, so most.charnock_viscosity = false turns the
  smooth-flow term off for Charnock/COARE but not for the others. Harmless while
  both defaults are true.
• Inputs.rst describes charnock_viscosity as "read only when
  roughness_type_sea = charnock"; it also gates coare3.0.
• The z₀ bounds are not documented in SurfaceLayer.rst. A flat 2.85 mm z₀ field
  in a tropical-cyclone run will look like a bug to anyone who plots it.
• adiabatic_donelan and adiabatic_wave_coupled still iterate on u* rather
  than z₀, and compare that velocity residual against tol_z, a length
  tolerance. Converges correctly; the other adiabatic variants were converted.
• Ri_b divides by θ while its numerator is a θ_v difference; t_star is defined
  from a dry-θ difference while ζ comes from a virtual-temperature Richardson
  number.
• Whitespace-only line at ERF_MOSTUtils.H:282.
• Source/PBL/ERF_ComputeDiffusivityMRF.cpp is whitespace-only in this diff and
  could be dropped to keep the PR focused.

---

Verification

ERF_MOSTStress.H and ERF_MOSTUtils.H compile clean under
g++ -std=c++20 -fsyntax-only -Wall -Wextra against the AMReX submodule headers
in both AMREX_USE_DOUBLE and AMREX_USE_FLOAT, with no warnings. The inner z₀
fixed point was checked numerically over z_ref ∈ {10, 20} m, depth ∈ {5, 10, 30,
100} m and U₁₀ ∈ {10 … 60} m/s: all cases converge, worst case 24 iterations.

This is a syntax and numerics check on the surface-layer headers, not a full
build, link, or regression run. The gold files in section 6 still need to be
exercised.

---

References

• Charnock, H. (1955), Wind stress on a water surface, QJRMS 81, 639–640.
• Smith, S. D. (1988), Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature, JGR 93, 15467–15472.
• Andreas, E. L. (1989), Thermal and size evolution of sea spray droplets, CRREL Rep. 89-11.
• Large, W. G., and S. Pond (1981), Open ocean momentum flux measurements in moderate to strong winds, JPO 11, 324–336.
• Beljaars, A. C. M. (1995), The parametrization of surface fluxes in large-scale models under free convection, QJRMS 121, 255–270.
• Fairall, C. W., et al. (2003), Bulk parameterization of air–sea fluxes: updates and verification for the COARE algorithm, J. Climate 16, 571–591.
• Grachev, A. A., C. W. Fairall, and E. F. Bradley (2000), Convective profile constants revisited, BLM 94, 495–515.
• Taylor, P. K., and M. J. Yelland (2001), The dependence of sea surface roughness on the height and steepness of the waves, JPO 31, 572–590.
• Donelan, M. A., et al. (2004), On the limiting aerodynamic roughness of the ocean in very strong winds, GRL 31, L18306, doi:10.1029/2004GL019460.
• Cheng, Y., and W. Brutsaert (2005), Flux-profile relationships for wind speed and temperature in the stable atmospheric boundary layer, BLM 114, 519–538.
• Davis, C., et al. (2008), Prediction of landfalling hurricanes with the Advanced Hurricane WRF model, MWR 136, 1990–2005, doi:10.1175/2007MWR2085.1.
• Warner, J. C., et al. (2010), Development of a coupled ocean–atmosphere–wave–sediment transport (COAWST) modeling system, Ocean Modelling 35, 230–244, doi:10.1016/j.ocemod.2010.07.010.
• Jiménez, P. A., et al. (2012), A revised scheme for the WRF surface layer formulation, MWR 140, 898–918, doi:10.1175/MWR-D-11-00056.1.
• Jiménez, P. A., and J. Dudhia (2018), On the need to modify the sea surface roughness formulation over shallow waters, JAMC 57, 1101–1110, doi:10.1175/JAMC-D-17-0137.1.

[AMLattanzi]

CPU nightly regression test failures:

MetGrid
WPS_Test_Terrain
WPS_Test_Terrain-OMP