pure subroutine drain_swept_flux_y(nx, ny, nz, areaT, uhh, hprev, aL, aR, a6, F)
!! Meridional analogue of drain_swept_flux_x.
integer, intent(in) :: nx, ny, nz
real(wp), intent(in) :: areaT(nx, ny)
real(wp), intent(in) :: uhh(nx, ny + 1, nz)
real(wp), intent(in) :: hprev(nx, ny, nz)
real(wp), intent(in) :: aL(nx, ny, nz), aR(nx, ny, nz), a6(nx, ny, nz)
real(wp), intent(inout) :: F(nx, ny + 1, nz)
integer :: i, j, k
real(wp) :: u, cfl, vol, conc
! Array-edge faces (1, ny+1) have no donor cell (cell 0 / ny+1 are
! off-array); they feed only re-wrapped ghosts, so zero them.
! Explicit do concurrent (not F(:,1,:)=0 array syntax -- offload).
do concurrent(k=1:nz, i=1:nx)
F(i, 1, k) = 0.0_wp
F(i, ny + 1, k) = 0.0_wp
end do
! Interior faces 2..ny: donors j-1 >= 1 (u>0) and j <= ny (u<0) in
! range. (Was j=1:ny+1, reading cell 0 / ny+1 -- latent OOB.)
do concurrent(k=1:nz, j=2:ny, i=1:nx) local(u, cfl, vol, conc)
u = uhh(i, j, k)
if (u > 0.0_wp) then
vol = max(areaT(i, j - 1)*hprev(i, j - 1, k), DRAIN_MIN_VOL)
cfl = min(u/vol, 1.0_wp)
conc = aR(i, j - 1, k) - 0.5_wp*cfl* &
((aR(i, j - 1, k) - aL(i, j - 1, k)) &
- a6(i, j - 1, k)*(1.0_wp - (2.0_wp/3.0_wp)*cfl))
F(i, j, k) = u*conc
else if (u < 0.0_wp) then
vol = max(areaT(i, j)*hprev(i, j, k), DRAIN_MIN_VOL)
cfl = min(-u/vol, 1.0_wp)
conc = aL(i, j, k) + 0.5_wp*cfl* &
((aR(i, j, k) - aL(i, j, k)) &
+ a6(i, j, k)*(1.0_wp - (2.0_wp/3.0_wp)*cfl))
F(i, j, k) = u*conc
else
F(i, j, k) = 0.0_wp
end if
end do
end subroutine drain_swept_flux_y