Fix cvmix_KPP StokesMOST: divide-by-zero and NaN injection under calm wind conditions

[patrickscholz]

Summary

Affects only cvmix_KPP with kpp_use_StokesMOST = .true. (the default). When wind speed is exactly zero, compute_stokes_velocities_MOM6style produces a division by zero and a 0/0 = NaN that permanently corrupts the Stokes drift and Langmuir turbulence in the affected column. The fix is a single-line lower bound on uv10.

---

The bug: a cascade of singularities at zero wind

The subroutine derives the peak wave frequency and wavenumber from the 10-m wind speed, then uses the wind direction to decompose the Stokes drift into vector components:

uv10    = sqrt(uwind**2 + vwind**2)   ! no floor — can be exactly 0.0

Tp      = 0.8 * uv10                  ! Tp = 0.0  (zero wave period)
fp      = 1.0 / Tp                    ! fp = Inf  (infinite frequency)
kp      = (2.0*pi*fp)**2 / g          ! kp = Inf  (infinite wavenumber → all waves at surface)

uS0_mag = Astokes * uv10**2           ! uS0_mag = 0.0
uS0     = uS0_mag * uwind / uv10      ! 0.0 / 0.0 = NaN
vS0     = uS0_mag * vwind / uv10      ! 0.0 / 0.0 = NaN

Every output of the subroutine (uS_top, vS_top, uS_mid, vS_mid, uS_mean, vS_mean) is then computed from uS0/vS0, so all six outputs become NaN. These NaN values feed directly into:

• cvmix_kpp_compute_StokesXi — the Stokes similarity parameter (Langmuir turbulence strength)
• compute_Efactor — the Langmuir enhancement factor for entrainment
• cvmix_coeffs_kpp — the final KPP diffusivities and viscosities

Once NaN enters the KPP mixing tensor it propagates irreversibly to tracer and momentum fields for that column — the model either crashes on the next timestep or silently produces garbage values.

When does zero wind occur?

• Doldrums / ITCZ: the equatorial calm belt where winds are climatologically weak and routinely pass through zero
• Idealised test cases: many standard configurations initialise with zero wind or prescribed surface fluxes without a wind field
• Temporal interpolation of forcing data: monthly-mean or 6-hourly forcing files are linearly interpolated between timesteps; any component that changes sign crosses through zero at each crossing
• Polar night / sea-ice edges: atmospheric forcing can be very small in winter polar regions

---

The fix

A single-line lower bound on uv10 before any ratio is computed (src/cvmix_driver/gen_modules_cvmix_kpp.F90):

! Before (wrong)
uv10 = sqrt(uwind**2 + vwind**2)

! After (correct)
uv10 = max(1.0e-3_WP, sqrt(uwind**2 + vwind**2))

1.0e-3 m/s (1 mm/s) represents physical background gustiness and is orders of magnitude below any meaningful wave-generating wind speed. It eliminates both the fp = 1/Tp singularity and the uwind/uv10 = 0/0 NaN while having zero practical effect on results.

---

Scope

• Only affects cvmix_KPP with kpp_use_StokesMOST = .true. (the default Langmuir turbulence path)
• No effect on runs without KPP or with kpp_use_StokesMOST = .false.
• No namelist changes required
• The fix is triggered only when sqrt(uwind^2 + vwind^2) < 1e-3 m/s; all other conditions are completely unaffected