Fix cvmix_KPP StokesMOST: OBL physical depth (metres) used as vertical level index

[patrickscholz]

Summary

Affects only cvmix_KPP with kpp_use_StokesMOST = .true. (the default). Two lines in the second-pass Stokes drift surface layer search were using kpp_obldepth (physical depth in metres) where kpp_nzobldepth (vertical level index) is required. This causes an out-of-bounds array access on every timestep in every ocean column where the OBL is deeper than the number of vertical levels.

---

Background: two companion OBL variables

cvmix_kpp_compute_OBL_depth produces two outputs that are always kept in sync:

Variable	Type	Meaning	Typical value
kpp_obldepth	real	OBL depth in metres	47.3 m
kpp_nzobldepth	real	OBL base as level index (+fraction)	5.7 (between levels 5 and 6)

These are fundamentally different quantities. kpp_obldepth is used for depth comparisons and physical diagnostics. kpp_nzobldepth is used whenever a level index is needed.

---

The bug

After the OBL depth is computed, KPP runs a second Stokes drift pass to find nzsfc — the level index of the surface layer base:

sldepth = kpp_surf_layer_ext * max(kpp_obldepth(node), kpp_minOBLdepth)  ! correct: metres
nzsfc   = kpp_obldepth(node)                ! WRONG: assigns ~47.3 m as a level index → nzsfc = 47
do nztmp = nun, int(kpp_obldepth(node))     ! WRONG: loop runs to level 47, array has ~5–10 levels in OBL
    if (-zbar_3d_n(nztmp+1,node) >= sldepth) then
        nzsfc = nztmp
        exit
    end if
end do

A typical FESOM2 production mesh has 47 vertical levels total. The OBL depth in the open ocean regularly reaches 50–500 m, meaning int(kpp_obldepth) produces values of 50–500. The loop then accesses zbar_3d_n(nztmp+1, node) for nztmp values far beyond the array's allocated size.

Why it doesn't always crash

Without compiler bounds checking, Fortran out-of-bounds reads return whatever memory happens to occupy that address. The loop's exit condition (-zbar >= sldepth) is eventually satisfied by chance on the garbage values, so nzsfc gets assigned a seemingly plausible but physically wrong level index. The model continues silently with a corrupted surface layer index feeding into the Stokes drift and Langmuir turbulence computation.

With bounds checking enabled (debug/development builds), this produces an immediate segfault.

---

The fix

Replace both erroneous uses of kpp_obldepth with kpp_nzobldepth (src/cvmix_driver/gen_modules_cvmix_kpp.F90):

! Before (wrong — physical depth in metres used as level index)
nzsfc = kpp_obldepth(node)
do nztmp = nun, int(kpp_obldepth(node))

! After (correct — level index companion used)
nzsfc = int(kpp_nzobldepth(node))
do nztmp = nun, int(kpp_nzobldepth(node))

The sldepth line immediately above is unchanged — it correctly uses kpp_obldepth for a physical depth comparison in metres.

---

Scope

• Only affects cvmix_KPP with kpp_use_StokesMOST = .true. (the default when using the StokesMOST Langmuir turbulence package)
• No effect on runs without KPP or with kpp_use_StokesMOST = .false.
• The corrupted nzsfc affects the Stokes drift profile and Langmuir enhancement factor in the second-pass surface layer computation — the primary path for Langmuir turbulence in the KPP boundary layer
• No namelist changes required