pure subroutine compute_channel_drag_rates(lambda_u, lambda_v, u_face, v_face, &
h_layer, wet_q, dyCu, dxCv, &
channel_drag, cdrag_side, &
nx_u, ny_u, nx_v, ny_v, nx, ny, nz)
!! Device kernel for the per-layer side-drag Rayleigh rate. See
!! `ocean_channel_drag_compute_tendencies` for the derivation.
!! `wet_q` is `(nx+1, ny+1)`; `dyCu` is `(nx+1, ny)` (u-face length
!! normal to the zonal flow); `dxCv` is `(nx, ny+1)`.
integer, intent(in) :: nx_u, ny_u, nx_v, ny_v, nx, ny, nz
logical, intent(in) :: channel_drag
real(wp), intent(in) :: cdrag_side
real(wp), intent(in) :: u_face(nx_u, ny_u, nz), v_face(nx_v, ny_v, nz)
real(wp), intent(in) :: h_layer(nx, ny, nz)
real(wp), intent(in) :: wet_q(nx + 1, ny + 1)
real(wp), intent(in) :: dyCu(nx + 1, ny), dxCv(nx, ny + 1)
real(wp), intent(out) :: lambda_u(nx_u, ny_u, nz), lambda_v(nx_v, ny_v, nz)
integer :: i, j, k
real(wp) :: f_s, f_n, f_blocked, speed, width, v_at_u, u_at_v
real(wp), parameter :: WIDTH_EPS = 1.0e-20_wp
! ---- Zero everywhere first (also the all-disabled short-circuit) ----
do concurrent(k=1:nz, j=1:ny_u, i=1:nx_u)
lambda_u(i, j, k) = 0.0_wp
end do
do concurrent(k=1:nz, j=1:ny_v, i=1:nx_v)
lambda_v(i, j, k) = 0.0_wp
end do
if (.not. channel_drag) return
if (cdrag_side <= 0.0_wp) return
! ---- u-faces: blocked by the two flanking corners (i,j),(i,j+1) ----
! Cross-stream cells at corner (i,j): T(i-1,j-1),T(i,j-1).
! Cross-stream cells at corner (i,j+1): T(i-1,j),T(i,j).
! Face range `j=2:ny` keeps the diagonal T reads in bounds: the
! south corner reads row `j-1` (>=1) and the north corner reads
! row `j` (<=ny); `wet_q(i,j+1)` reaches `j+1<=ny+1`, in bounds.
! Row `j=1` is a ghost row with no prognostic velocity, so leaving
! its rate at 0 is exact.
do concurrent(k=1:nz, j=2:ny, i=2:nx) &
local(f_s, f_n, f_blocked, speed, width, v_at_u)
! South corner (i,j): land OR layer vanished in T(i-1,j-1),T(i,j-1).
f_s = 1.0_wp - wet_q(i, j)
if (min(h_layer(i - 1, j - 1, k), h_layer(i, j - 1, k)) < SIDE_H_VANISH) then
f_s = 1.0_wp
end if
! North corner (i,j+1): land OR layer vanished in T(i-1,j),T(i,j).
f_n = 1.0_wp - wet_q(i, j + 1)
if (min(h_layer(i - 1, j, k), h_layer(i, j, k)) < SIDE_H_VANISH) then
f_n = 1.0_wp
end if
f_blocked = 0.5_wp*f_s + 0.5_wp*f_n
v_at_u = 0.25_wp*( &
v_face(i - 1, j, k) + v_face(i, j, k) + &
v_face(i - 1, j + 1, k) + v_face(i, j + 1, k))
speed = sqrt(u_face(i, j, k)*u_face(i, j, k) + v_at_u*v_at_u)
width = max(dyCu(i, j), WIDTH_EPS)
lambda_u(i, j, k) = cdrag_side*speed*f_blocked/width
end do
! ---- v-faces: blocked by the two flanking corners (i,j),(i+1,j) ----
! Cross-stream cells at corner (i,j): T(i-1,j-1),T(i-1,j).
! Cross-stream cells at corner (i+1,j): T(i,j-1),T(i,j).
! Face range `i=2:nx` keeps the diagonal T reads in bounds: the
! west corner reads column `i-1` (>=1) and the east corner reads
! column `i` (<=nx); `wet_q(i+1,j)` reaches `i+1<=nx+1`, in bounds.
! Column `i=1` is a ghost column with no prognostic velocity.
do concurrent(k=1:nz, j=2:ny, i=2:nx) &
local(f_s, f_n, f_blocked, speed, width, u_at_v)
! West corner (i,j): land OR layer vanished in T(i-1,j-1),T(i-1,j).
f_s = 1.0_wp - wet_q(i, j)
if (min(h_layer(i - 1, j - 1, k), h_layer(i - 1, j, k)) < SIDE_H_VANISH) then
f_s = 1.0_wp
end if
! East corner (i+1,j): land OR layer vanished in T(i,j-1),T(i,j).
f_n = 1.0_wp - wet_q(i + 1, j)
if (min(h_layer(i, j - 1, k), h_layer(i, j, k)) < SIDE_H_VANISH) then
f_n = 1.0_wp
end if
f_blocked = 0.5_wp*f_s + 0.5_wp*f_n
u_at_v = 0.25_wp*( &
u_face(i, j - 1, k) + u_face(i + 1, j - 1, k) + &
u_face(i, j, k) + u_face(i + 1, j, k))
speed = sqrt(v_face(i, j, k)*v_face(i, j, k) + u_at_v*u_at_v)
width = max(dxCv(i, j), WIDTH_EPS)
lambda_v(i, j, k) = cdrag_side*speed*f_blocked/width
end do
end subroutine compute_channel_drag_rates