region_power Subroutine

   subroutine region_power(grid, ms, bt_work, F_u_3d, F_v_3d, j_lo, j_hi, p_mean)
      !! Cell-centred BT power per region: P = ⟨u_bt·F_u + v_bt·F_v⟩.
      !! `F_u_3d` is per-layer at u-faces, shape `(nx+1, ny, nz)`.
      !! Depth-averages with the centred face thickness as weight to
      !! get the BT-mode contribution, then dots with `bt_ubt/bt_vbt`
      !! at the same face, then averages east+west (north+south) into
      !! the cell.
      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) :: F_u_3d(:, :, :), F_v_3d(:, :, :)
      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, 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
            ! Depth-mean each tendency at the four cell faces.  Weight
            ! by the centred face thickness so the BT mode is what we
            ! actually project onto.
            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 + F_u_3d(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 + F_u_3d(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 + F_v_3d(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 + F_v_3d(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
            ! Power per face: u_bt·F.  Average into the cell.
            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