vmix_assemble_clip_henyey_impl Subroutine

   pure subroutine vmix_assemble_clip_henyey_impl(nx, ny, nzp1, kv, kt, ks, &
                                                  kv_bg, kt_bg, ks_bg, kv_max, kd_max, &
                                                  n0_2omega, henyey_max_lat, kd_min, &
                                                  geolat)
      !! `vmix_assemble_clip_impl` with the Henyey latitude factor
      !! (`henyey_lat_factor_impl`) scaling the two SCALAR TRACER background
      !! floors, each then floored at the minimum diffusivity:
      !!
      !!   kt floor = max(kd_min, kt_bg · L(φ))
      !!   ks floor = max(kd_min, ks_bg · L(φ))
      !!
      !! matching the reference `Kd_sfc = max(Kd_min, Kd · L)`.  Without the
      !! `kd_min` floor the background would collapse toward zero at the
      !! equator (`L(0°) = 0` exactly) and at the poleward clamp.
      !!
      !! The momentum floor stays the plain scalar `kv_bg` — see the
      !! `bkgnd_henyey` docstring for why (rdb's `kv_bg` is an
      !! independent momentum floor, not `prandtl · kt_bg`).
      !!
      !! Bryan-Lewis composition is impossible here by construction: this
      !! routine takes no depth profile, and `vmix_assemble` reaches it only
      !! on the `bkgnd_profile == .false.` arm.
      !!
      !! Interior interfaces k = 2..nzp1-1.  The latitude factor is a column
      !! constant, so it is computed once per (i,j) and reused down the
      !! column.  No optional arguments: the Henyey parameter group is
      !! supplied in full or this routine is not the one you call.
      !!
      !! Explicit-shape args so the do concurrent stays descriptor-walk free.
      integer, intent(in) :: nx, ny, nzp1
      real(wp), intent(inout) :: kv(nx, ny, nzp1)
      real(wp), intent(inout) :: kt(nx, ny, nzp1)
      real(wp), intent(inout) :: ks(nx, ny, nzp1)
      real(wp), intent(in) :: kv_bg, kt_bg, ks_bg, kv_max, kd_max
      real(wp), intent(in) :: n0_2omega, henyey_max_lat
         !! `N0/(2Ω)` ratio (nondim) and the poleward cutoff latitude (degN).
      real(wp), intent(in) :: kd_min
         !! Minimum background tracer diffusivity (m^2/s), ALREADY resolved
         !! through `vmix_resolve_kd_min` — this kernel never sees the
         !! negative sentinel.
      real(wp), intent(in) :: geolat(nx, ny)
         !! T-point latitude (degN) — `ocean_metrics_t%geolatT`, the T
         !! stagger.  The corner field `geolatBu` is `(nx+1, ny+1)` and is
         !! NOT a substitute.

      integer :: i, j, k
      real(wp) :: kt_floor, ks_floor

      do concurrent(j=1:ny, i=1:nx) local(kt_floor, ks_floor, k)
         ! ONE write site per local (a `local()` variable reassigned across
         ! a branch is the recorded gfortran-15.1 corruption shape).
         kt_floor = max(kd_min, kt_bg* &
                        henyey_lat_factor_impl(geolat(i, j), n0_2omega, henyey_max_lat))
         ks_floor = max(kd_min, ks_bg* &
                        henyey_lat_factor_impl(geolat(i, j), n0_2omega, henyey_max_lat))
         do k = 2, nzp1 - 1
            kv(i, j, k) = min(max(kv(i, j, k), kv_bg), kv_max)
            kt(i, j, k) = min(max(kt(i, j, k), kt_floor), kd_max)
            ks(i, j, k) = min(max(ks(i, j, k), ks_floor), kd_max)
         end do
      end do
   end subroutine vmix_assemble_clip_henyey_impl