subroutine region_power_drag_implicit(grid, ms, bt_work, rate_u, rate_v, &
j_lo, j_hi, p_mean)
!! Bottom drag in implicit mode: the slow tendency for u at a
!! face is `−rate_u(i,j,k)·u_face_layer(i,j,k)` (1/s × m/s
!! → m/s²). We construct that on the fly, depth-mean by
!! face thickness, dot with `bt_ubt` (and v counterpart),
!! sum the cell-centred result over the region. Sign comes
!! out negative: drag removes BT-mode energy.
type(hgrid_t), intent(in) :: grid
type(multilayer_state_t), intent(in) :: ms
type(barotropic_workstate_t), intent(in) :: bt_work
real(wp), intent(in) :: rate_u(:, :, :), rate_v(:, :, :)
integer, intent(in) :: j_lo, j_hi
real(wp), intent(out) :: p_mean
integer :: i, j, k, nz, ip_lo, ip_hi, n
real(wp) :: h_face
real(wp) :: fu_W, fu_E, fv_S, fv_N
real(wp) :: w_W, w_E, w_S, w_N
real(wp) :: power_sum, p_cell
ip_lo = grid%nghost + 1
ip_hi = grid%nghost + grid%nx_phys
nz = ms%nz_ml
power_sum = 0.0_wp
n = 0
do j = j_lo, j_hi
do i = ip_lo, ip_hi
fu_W = 0.0_wp
w_W = 0.0_wp
fu_E = 0.0_wp
w_E = 0.0_wp
fv_S = 0.0_wp
w_S = 0.0_wp
fv_N = 0.0_wp
w_N = 0.0_wp
do k = 1, nz
h_face = 0.5_wp*(ms%h_layer(i - 1, j, k) + ms%h_layer(i, j, k))
fu_W = fu_W + (-rate_u(i, j, k)*ms%u_face_x_layer(i, j, k))*h_face
w_W = w_W + h_face
h_face = 0.5_wp*(ms%h_layer(i, j, k) + ms%h_layer(i + 1, j, k))
fu_E = fu_E + (-rate_u(i + 1, j, k)*ms%u_face_x_layer(i + 1, j, k))*h_face
w_E = w_E + h_face
h_face = 0.5_wp*(ms%h_layer(i, j - 1, k) + ms%h_layer(i, j, k))
fv_S = fv_S + (-rate_v(i, j, k)*ms%v_face_y_layer(i, j, k))*h_face
w_S = w_S + h_face
h_face = 0.5_wp*(ms%h_layer(i, j, k) + ms%h_layer(i, j + 1, k))
fv_N = fv_N + (-rate_v(i, j + 1, k)*ms%v_face_y_layer(i, j + 1, k))*h_face
w_N = w_N + h_face
end do
if (w_W > 0.0_wp) fu_W = fu_W/w_W
if (w_E > 0.0_wp) fu_E = fu_E/w_E
if (w_S > 0.0_wp) fv_S = fv_S/w_S
if (w_N > 0.0_wp) fv_N = fv_N/w_N
p_cell = 0.5_wp*(bt_work%bt_ubt(i, j)*fu_W + bt_work%bt_ubt(i + 1, j)*fu_E) + &
0.5_wp*(bt_work%bt_vbt(i, j)*fv_S + bt_work%bt_vbt(i, j + 1)*fv_N)
power_sum = power_sum + p_cell
n = n + 1
end do
end do
if (n > 0) then
p_mean = power_sum/real(n, wp)
else
p_mean = 0.0_wp
end if
end subroutine region_power_drag_implicit