pure subroutine gm_column_y(nx, ny, nz, i_smax2, i4dt, use_open, &
dx_cv, areaT, h_layer, bathy, slope_y, khth_v, open_v, vhD)
!! v-face GM column recurrence — mirror of `gm_column_x` with the
!! v-stagger. Interior v-face (j=2..ny) pairs columns js=j-1 (south)
!! and j (north). See `gm_column_x` for the open-column (`use_open`)
!! rule and the MOM6 `nk_linear` divergence note.
integer, intent(in) :: nx, ny, nz
real(wp), intent(in) :: i_smax2, i4dt
real(wp), intent(in) :: dx_cv(nx, ny + 1)
real(wp), intent(in) :: areaT(nx, ny)
real(wp), intent(in) :: h_layer(nx, ny, nz)
real(wp), intent(in) :: bathy(nx, ny)
real(wp), intent(in) :: slope_y(nx, ny + 1, nz + 1)
real(wp), intent(in) :: khth_v(nx, ny + 1)
logical, intent(in) :: use_open
real(wp), intent(in) :: open_v(nx, ny + 1, nz)
real(wp), intent(inout) :: vhD(nx, ny + 1, nz)
integer :: i, j, k, js, ka, kb, ktop
real(wp) :: havS(NZ_STACK_MAX), havN(NZ_STACK_MAX)
real(wp) :: rsumS(NZ_STACK_MAX + 1), rsumN(NZ_STACK_MAX + 1)
real(wp) :: eS(NZ_STACK_MAX + 1), eN(NZ_STACK_MAX + 1)
logical :: ok(NZ_STACK_MAX)
real(wp) :: vhtot, slope, s2r, sfn_unlim, sfn_safe, sfn_est, sfn_in_h
real(wp) :: h_frac_d, vhd_k, kh
do concurrent(j=2:ny, i=1:nx) &
local(k, js, ka, kb, ktop, havS, havN, rsumS, rsumN, eS, eN, ok, vhtot, slope, &
s2r, sfn_unlim, sfn_safe, sfn_est, sfn_in_h, h_frac_d, vhd_k, kh)
js = j - 1
kh = khth_v(i, j)
eS(1) = -bathy(i, js)
eN(1) = -bathy(i, j)
do k = 1, nz
eS(k + 1) = eS(k) + h_layer(i, js, k)
eN(k + 1) = eN(k) + h_layer(i, j, k)
end do
do k = 1, nz
ok(k) = .true.
end do
if (use_open) then
do k = 1, nz
ok(k) = open_v(i, j, k) > 0.5_wp .and. &
rdb_vl_is_live(h_layer(i, js, k)) .and. &
rdb_vl_is_live(h_layer(i, j, k))
end do
end if
ktop = 0
do k = nz, 1, -1
if (ok(k)) then
ktop = k
exit
end if
end do
rsumS(nz + 1) = 0.0_wp
rsumN(nz + 1) = 0.0_wp
do k = nz, 1, -1
if (ok(k)) then
havS(k) = max(i4dt*areaT(i, js)*(h_layer(i, js, k) - H_VANISHED), 0.0_wp)
havN(k) = max(i4dt*areaT(i, j)*(h_layer(i, j, k) - H_VANISHED), 0.0_wp)
else
havS(k) = 0.0_wp
havN(k) = 0.0_wp
end if
rsumS(k) = rsumS(k + 1) + havS(k)
rsumN(k) = rsumN(k + 1) + havN(k)
end do
vhtot = 0.0_wp
vhD(i, j, nz) = 0.0_wp
do k = 2, nz
ka = k
kb = k - 1
if (.not. ok(kb) .or. kb >= ktop) then
vhD(i, j, kb) = 0.0_wp
cycle
end if
slope = slope_y(i, j, k)
if (.not. ieee_is_finite(slope)) slope = 0.0_wp
s2r = slope*slope*i_smax2
sfn_unlim = gm_block_below_bed(-(kh*dx_cv(i, j))*slope, &
eS(k), eS(k - 1), eN(1), eN(k), eN(k - 1), eS(1))
if (vhtot <= 0.0_wp) then
h_frac_d = gm_h_frac(havS(kb), rsumS(kb))
else
h_frac_d = gm_h_frac(havN(kb), rsumN(kb))
end if
sfn_safe = vhtot*(1.0_wp - h_frac_d)
sfn_est = (sfn_unlim + s2r*sfn_safe)/(1.0_wp + s2r)
sfn_in_h = min(max(sfn_est, -rsumS(ka)), rsumN(ka))
vhd_k = max(min(sfn_in_h - vhtot, havS(kb)), -havN(kb))
vhD(i, j, kb) = vhd_k
vhtot = vhtot + vhd_k
end do
if (ktop > 0) vhD(i, j, ktop) = -vhtot
end do
end subroutine gm_column_y