henyey_lat_factor_impl Function

   pure function henyey_lat_factor_impl(lat_deg, n0_2omega, max_lat) result(fac)
      !$acc routine seq
      !! Henyey, Wright & Flatte (1986) JGR 91:8487 latitude dependence of
      !! the internal-wave-driven mixing rate, in the SIMPLIFIED constant-`N0`
      !! form of Harrison & Hallberg (2008) JPO 38:1894 — the in-situ column
      !! stratification is replaced by a fixed reference `N0`, so the factor
      !! collapses to a pure function of latitude:
      !!
      !!   L(φ) = |sin φ| · acosh(N0_2Ω / max(|sin φ|, HENYEY_MIN_SINLAT))
      !!          / [ sin 30° · acosh(N0_2Ω / sin 30°) ]
      !!
      !! `N0_2Ω = N0/(2Ω)`.  `L(30°) = 1` to round-off by construction (the
      !! denominator IS the numerator at 30°), which is what makes 30°
      !! useless as a test latitude for "was the factor applied at all" —
      !! a test that cannot distinguish `×L` from `×1` there.  Use 45°/90°.
      !!
      !! Equator singularity: `|sin φ|` is floored to `HENYEY_MIN_SINLAT`
      !! ONLY inside the ratio that feeds `acosh` (that is the 1/0 guard);
      !! the OUTER multiplication keeps the TRUE (unfloored) `|sin φ|`, so
      !! `L → 0` smoothly at the equator instead of blowing up, and
      !! `L(0°) = 0` EXACTLY.
      !!
      !! Poleward of `max_lat` (degN, compared against `|lat_deg|` so BOTH
      !! hemispheres clamp) `|sin φ|` is reset to `HENYEY_MIN_SINLAT`
      !! everywhere in the expression, collapsing `L` to a tiny positive
      !! floor (~1.2e-9 at the default `n0_2omega`) — NOT to the exact zero
      !! the equator produces.  Inert at the default `max_lat = 95` (> 90,
      !! never trips for a real latitude).
      !!
      !! This function returns the RAW factor `L(phi)` — it applies no
      !! minimum-diffusivity floor, so `L(0°)` really is exactly 0 and the
      !! poleward clamp really does return ~1.2e-9.  The floor lives one
      !! level up, in `vmix_assemble_clip_henyey_impl`, as
      !! `max(kd_min, kt_bg·L(phi))` — matching the reference code, which
      !! wraps the scaled diffusivity in `max(Kd_min, Kd·L)` with `Kd_min`
      !! defaulting to `0.01·Kd` (`HENYEY_KD_MIN_FRAC`).  Keeping the floor
      !! out of the factor is what lets the unit tests assert the closed
      !! form of `L` directly against an `acosh` oracle.
      !!
      !! `pure` + `!$acc routine seq`: called per column from the
      !! `do concurrent` in `vmix_assemble_clip_henyey_impl`, and directly
      !! from the unit tests as the analytical oracle's subject.
      real(wp), intent(in) :: lat_deg
         !! T-point latitude, degrees north (may be negative).
      real(wp), intent(in) :: n0_2omega
         !! `N0/(2Ω)` reference-stratification ratio (nondim, ≥ 1).
      real(wp), intent(in) :: max_lat
         !! Poleward cutoff latitude (degN, compared against `|lat_deg|`).
      real(wp) :: fac

      real(wp) :: abs_sinlat

      ! ONE write site per local (a `local()` variable reassigned across a
      ! branch is the recorded gfortran-15.1 corruption shape) — `merge`
      ! evaluates both arms, and both are finite for every input: |sin| is
      ! bounded by 1 and n0_2omega ≥ 1, so the acosh argument is always ≥ 1.
      abs_sinlat = merge(HENYEY_MIN_SINLAT, abs(sin(lat_deg*DEG2RAD)), &
                         abs(lat_deg) > max_lat)
      ! sin(30 deg) = 0.5 exactly; the denominator normalises L(30 deg) = 1.
      fac = abs_sinlat*acosh(n0_2omega/max(HENYEY_MIN_SINLAT, abs_sinlat)) &
            /(0.5_wp*acosh(2.0_wp*n0_2omega))
   end function henyey_lat_factor_impl