pure subroutine varmix_sn_v(nx, ny, nz, do_visbeck, s2max, h_layer, &
slope_x, slope_y, n2_v, sn_v)
!! Thickness-weighted Eady growth rate at v-faces (mirror of
!! `varmix_sn_u`). Interior v-face (j=2..ny) pairs js=j-1 (south) + j.
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_v(nx, ny + 1, nz + 1)
real(wp), intent(out) :: sn_v(nx, ny + 1)
integer :: i, j, k, js, im, ip
real(wp) :: hgeom, hdn, hup, s2, n2, sn_acc, h_acc
real(wp) :: wsw, wse, wnw, wne, sx2, denom
do concurrent(j=1:ny + 1, i=1:nx) &
local(k, js, im, ip, hgeom, hdn, hup, s2, n2, sn_acc, h_acc, &
wsw, wse, wnw, wne, sx2, denom)
sn_v(i, j) = 0.0_wp
if (do_visbeck .and. j >= 2 .and. j <= ny) then
js = j - 1
im = max(1, i - 1) ! west cell-column (array-edge clamp)
ip = min(nx, i + 1) ! east cell-column (array-edge clamp)
sn_acc = 0.0_wp
h_acc = 0.0_wp
do k = 2, nz
hdn = sqrt(max(h_layer(i, js, k)*h_layer(i, j, k), 0.0_wp))
hup = sqrt(max(h_layer(i, js, k - 1)*h_layer(i, j, k - 1), 0.0_wp))
hgeom = sqrt(hdn*hup)
! 4 corner slope_x values around the v-face, weighted by the
! MOM6 h4_u product co-located with each corner's slope_x
! u-point: the 4 thicknesses straddling that u-face (the two
! cell-cols it separates) at the two layers (k, k-1) the
! interface separates. Uniform h ⇒ bit-identical.
wse = (h_layer(i, js, k)*h_layer(ip, js, k)) &
*(h_layer(i, js, k - 1)*h_layer(ip, js, k - 1))
wnw = (h_layer(im, j, k)*h_layer(i, j, k)) &
*(h_layer(im, j, k - 1)*h_layer(i, j, k - 1))
wne = (h_layer(i, j, k)*h_layer(ip, j, k)) &
*(h_layer(i, j, k - 1)*h_layer(ip, j, k - 1))
wsw = (h_layer(im, js, k)*h_layer(i, js, k)) &
*(h_layer(im, js, k - 1)*h_layer(i, js, k - 1))
denom = ((wse + wnw) + (wne + wsw)) + H_SUBROUNDOFF4
sx2 = (((wse*slope_x(i + 1, js, k)*slope_x(i + 1, js, k)) + &
(wnw*slope_x(i, j, k)*slope_x(i, j, k))) + &
((wne*slope_x(i + 1, j, k)*slope_x(i + 1, j, k)) + &
(wsw*slope_x(i, js, k)*slope_x(i, js, k))))/denom
s2 = slope_y(i, j, k)*slope_y(i, j, k) + sx2
if (s2max > 0.0_wp) s2 = s2*s2max/(s2 + s2max)
n2 = max(n2_v(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_v(i, j) = sn_acc/h_acc
end if
end do
end subroutine varmix_sn_v