varmix_sn_u Subroutine

   pure subroutine varmix_sn_u(nx, ny, nz, do_visbeck, s2max, h_layer, &
                               slope_x, slope_y, n2_u, sn_u)
      !! Thickness-weighted Eady growth rate at u-faces (own component).
      !! Interior u-face (i=2..nx) pairs centre columns iw=i-1 (west) and i
      !! (east).  Interior interfaces K=2..nz; `H_geom = sqrt(sqrt(h_iw,k *
      !! h_i,k) * sqrt(h_iw,k-1 * h_i,k-1))`.  S2 = slope_x^2 + the four
      !! corner slope_y^2 h-weighted to the u-face; S2 optionally limited.
      !! `SN_u = sum sqrt(S2*N2)*H_geom / sum H_geom`.
      integer, intent(in) :: nx, ny, nz
      logical, intent(in) :: do_visbeck
      real(wp), intent(in) :: s2max
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: slope_x(nx + 1, ny, nz + 1)
      real(wp), intent(in) :: slope_y(nx, ny + 1, nz + 1)
      real(wp), intent(in) :: n2_u(nx + 1, ny, nz + 1)
      real(wp), intent(out) :: sn_u(nx + 1, ny)

      integer :: i, j, k, iw, jm, jp
      real(wp) :: hgeom, hdn, hup, s2, n2, sn_acc, h_acc
      real(wp) :: wsw, wse, wnw, wne, sy2, denom

      do concurrent(j=1:ny, i=1:nx + 1) &
         local(k, iw, jm, jp, hgeom, hdn, hup, s2, n2, sn_acc, h_acc, &
               wsw, wse, wnw, wne, sy2, denom)
         sn_u(i, j) = 0.0_wp
         if (do_visbeck .and. i >= 2 .and. i <= nx) then
            iw = i - 1
            jm = max(1, j - 1)        ! south cell-row (array-edge clamp)
            jp = min(ny, j + 1)       ! north cell-row (array-edge clamp)
            sn_acc = 0.0_wp
            h_acc = 0.0_wp
            do k = 2, nz       ! interior interface index
               hdn = sqrt(max(h_layer(iw, j, k)*h_layer(i, j, k), 0.0_wp))
               hup = sqrt(max(h_layer(iw, j, k - 1)*h_layer(i, j, k - 1), 0.0_wp))
               hgeom = sqrt(hdn*hup)
               ! 4 corner slope_y values around the u-face, weighted by the
               ! MOM6 h4_v product co-located with each corner's slope_y
               ! v-point: the 4 thicknesses straddling that v-face (the two
               ! cell-rows it separates) at the two layers (k, k-1) the
               ! interface separates.  Under uniform h all four are equal ⇒
               ! bit-identical to a single-thickness weight.
               wnw = (h_layer(iw, j, k)*h_layer(iw, jp, k)) &
                     *(h_layer(iw, j, k - 1)*h_layer(iw, jp, k - 1))
               wne = (h_layer(i, j, k)*h_layer(i, jp, k)) &
                     *(h_layer(i, j, k - 1)*h_layer(i, jp, k - 1))
               wsw = (h_layer(iw, jm, k)*h_layer(iw, j, k)) &
                     *(h_layer(iw, jm, k - 1)*h_layer(iw, j, k - 1))
               wse = (h_layer(i, jm, k)*h_layer(i, j, k)) &
                     *(h_layer(i, jm, k - 1)*h_layer(i, j, k - 1))
               denom = ((wse + wnw) + (wne + wsw)) + H_SUBROUNDOFF4
               sy2 = (((wnw*slope_y(iw, j + 1, k)*slope_y(iw, j + 1, k)) + &
                       (wse*slope_y(i, j, k)*slope_y(i, j, k))) + &
                      ((wne*slope_y(i, j + 1, k)*slope_y(i, j + 1, k)) + &
                       (wsw*slope_y(iw, j, k)*slope_y(iw, j, k))))/denom
               s2 = slope_x(i, j, k)*slope_x(i, j, k) + sy2
               if (s2max > 0.0_wp) s2 = s2*s2max/(s2 + s2max)
               n2 = max(n2_u(i, j, k), 0.0_wp)
               sn_acc = sn_acc + sqrt(s2*n2)*hgeom
               h_acc = h_acc + hgeom
            end do
            if (h_acc > 0.0_wp) sn_u(i, j) = sn_acc/h_acc
         end if
      end do
   end subroutine varmix_sn_u