pure subroutine ocean_slopes_pass_y(nx, ny, nz, eos, rho0, min_dz, &
h_layer, t_fill, s_fill, e_int, &
idyCv, wet_v, slope_y, n2_v)
!! v-face slope + N² pass — mirror of `pass_x` with v-staggering.
!! The v-face at (i,j) sits between cells (i,j-1) and (i,j); pairs
!! columns `js=j-1` (south) and `j` (north), loop `j=2:ny`.
integer, intent(in) :: nx, ny, nz
type(eos_t), intent(in) :: eos
real(wp), intent(in) :: rho0, min_dz
real(wp), intent(in) :: h_layer(nx, ny, nz)
real(wp), intent(in) :: t_fill(nx, ny, nz)
real(wp), intent(in) :: s_fill(nx, ny, nz)
real(wp), intent(in) :: e_int(nx, ny, nz + 1)
real(wp), intent(in) :: idyCv(nx, ny + 1)
real(wp), intent(in) :: wet_v(nx, ny + 1)
real(wp), intent(out) :: slope_y(nx, ny + 1, nz + 1)
real(wp), intent(out) :: n2_v(nx, ny + 1, nz + 1)
integer :: i, j, k, js, ka, kb
real(wp) :: pres_v, t_v, s_v, rho_v, dsv_dt, dsv_ds, drdt, drds
real(wp) :: drdjA, drdjB, drdkL, drdkR
real(wp) :: hg2A, hg2B, hg2L, hg2R, haA, haB, haL, haR
real(wp) :: dzaL, dzaR, wtA, wtB, wtL, wtR
real(wp) :: drdy, drdz, mag2, slope, presS, presN
real(wp) :: g_rho0, mask
g_rho0 = GRAVITY/rho0
do concurrent(i=1:nx, j=1:ny + 1)
slope_y(i, j, 1) = 0.0_wp
slope_y(i, j, nz + 1) = 0.0_wp
n2_v(i, j, 1) = 0.0_wp
n2_v(i, j, nz + 1) = 0.0_wp
end do
do concurrent(k=1:nz + 1, i=1:nx)
slope_y(i, 1, k) = 0.0_wp
slope_y(i, ny + 1, k) = 0.0_wp
n2_v(i, 1, k) = 0.0_wp
n2_v(i, ny + 1, k) = 0.0_wp
end do
do concurrent(k=2:nz, j=2:ny, i=1:nx) &
local(js, ka, kb, pres_v, t_v, s_v, rho_v, dsv_dt, dsv_ds, &
drdt, drds, drdjA, drdjB, drdkL, drdkR, &
hg2A, hg2B, hg2L, hg2R, haA, haB, haL, haR, &
dzaL, dzaR, wtA, wtB, wtL, wtR, drdy, drdz, &
mag2, slope, presS, presN, mask)
js = j - 1
ka = k
kb = k - 1
presS = pressure_above_x(nx, ny, nz, h_layer, i, js, ka, rho0)
presN = pressure_above_x(nx, ny, nz, h_layer, i, j, ka, rho0)
pres_v = 0.5_wp*(presS + presN)
t_v = 0.25_wp*((t_fill(i, js, ka) + t_fill(i, j, ka)) + &
(t_fill(i, js, kb) + t_fill(i, j, kb)))
s_v = 0.25_wp*((s_fill(i, js, ka) + s_fill(i, j, ka)) + &
(s_fill(i, js, kb) + s_fill(i, j, kb)))
call eos_density_specvol_derivs(eos, t_v, s_v, pres_v, rho_v, dsv_dt, dsv_ds)
drdt = -(rho_v*rho_v)*dsv_dt
drds = -(rho_v*rho_v)*dsv_ds
drdjA = drdt*(t_fill(i, j, ka) - t_fill(i, js, ka)) + &
drds*(s_fill(i, j, ka) - s_fill(i, js, ka))
drdjB = drdt*(t_fill(i, j, kb) - t_fill(i, js, kb)) + &
drds*(s_fill(i, j, kb) - s_fill(i, js, kb))
drdkL = drdt*(t_fill(i, js, kb) - t_fill(i, js, ka)) + &
drds*(s_fill(i, js, kb) - s_fill(i, js, ka))
drdkR = drdt*(t_fill(i, j, kb) - t_fill(i, j, ka)) + &
drds*(s_fill(i, j, kb) - s_fill(i, j, ka))
hg2A = h_layer(i, js, ka)*h_layer(i, j, ka) + H_DIV_EPS*H_DIV_EPS
hg2B = h_layer(i, js, kb)*h_layer(i, j, kb) + H_DIV_EPS*H_DIV_EPS
hg2L = h_layer(i, js, ka)*h_layer(i, js, kb) + H_DIV_EPS*H_DIV_EPS
hg2R = h_layer(i, j, ka)*h_layer(i, j, kb) + H_DIV_EPS*H_DIV_EPS
haA = 0.5_wp*(h_layer(i, js, ka) + h_layer(i, j, ka)) + H_DIV_EPS
haB = 0.5_wp*(h_layer(i, js, kb) + h_layer(i, j, kb)) + H_DIV_EPS
haL = 0.5_wp*(h_layer(i, js, ka) + h_layer(i, js, kb)) + H_DIV_EPS
haR = 0.5_wp*(h_layer(i, j, ka) + h_layer(i, j, kb)) + H_DIV_EPS
dzaL = max(haL, min_dz)
dzaR = max(haR, min_dz)
wtA = hg2A*haB
wtB = hg2B*haA
wtL = hg2L*(haR*dzaR)
wtR = hg2R*(haL*dzaL)
drdz = ((wtL*drdkL) + (wtR*drdkR))/((dzaL*wtL) + (dzaR*wtR))
drdy = ((wtA*drdjA + wtB*drdjB)/(wtA + wtB) - &
drdz*(e_int(i, js, k) - e_int(i, j, k)))*idyCv(i, j)
mag2 = drdy*drdy + drdz*drdz
if (mag2 > 0.0_wp) then
slope = drdy/sqrt(mag2)
else
slope = 0.0_wp
end if
mask = wet_v(i, j)
slope_y(i, j, k) = slope*mask
n2_v(i, j, k) = g_rho0*drdz*mask
end do
end subroutine ocean_slopes_pass_y