compute_bt_rem Subroutine

   pure subroutine compute_bt_rem(grid, bt_work, ms, metrics, r_linear, hbbl, dt_inner)
      !! Per-face multiplicative damping factor for the BT-substep velocity
      !! update (linear-drag branch):
      !!     bt_rem_face = Htot_face / (Htot_face + r·hbbl·dt_inner)
      !! applied as ubt_new = bt_rem_u·(ubt_old + dt_inner·forces) each inner
      !! step. When the `bt_substep_drag` knob is off this must NOT be called and
      !! the workspace stays at 1 (no-op, bit-identical).
      type(hgrid_t), intent(in) :: grid
      type(barotropic_workstate_t), intent(inout) :: bt_work
      type(multilayer_state_t), intent(in) :: ms
      type(ocean_metrics_t), intent(in) :: metrics
         !! REQUIRED (see `derive_bt_from_layers`).  Under
         !! `&vcoord_nml zfixed_closed_faces` `Htot_face` is the OPEN-column
         !! face depth `Σ_k h_face·open` — the column `ubt` is the mean of
         !! and the one the fast loop transports on — so the bed stress
         !! `r·hbbl·ubt` is spread over the depth that actually carries the
         !! barotropic momentum.  The full-column depth would under-damp a
         !! partially closed face by its closed fraction.  `.false.` (the
         !! default) ⇒ the original loops, textually unchanged; `open_*`
         !! is never named.
      real(wp), intent(in) :: r_linear   !! Linear drag rate at the bed (1/s)
      real(wp), intent(in) :: hbbl       !! BBL thickness over which drag acts (m)
      real(wp), intent(in) :: dt_inner   !! BT-substep dt (s)

      integer :: i, j, k, nx, ny, nz
      real(wp) :: htot_face, drag_dt

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      drag_dt = r_linear*hbbl*dt_inner   ! product is in metres

      if (metrics%use_closed_faces) then
         call bt_rem_open_impl(nx, ny, nz, ms%h_layer, metrics%open_u, &
                               metrics%open_v, drag_dt, &
                               bt_work%bt_rem_u, bt_work%bt_rem_v)
         return
      end if

      do concurrent(j=1:ny, i=2:nx) local(k, htot_face)
         htot_face = 0.0_wp
         do k = 1, nz
            htot_face = htot_face + 0.5_wp*(ms%h_layer(i - 1, j, k) + ms%h_layer(i, j, k))
         end do
         if (htot_face > 0.0_wp) then
            bt_work%bt_rem_u(i, j) = htot_face/(htot_face + drag_dt)
         else
            bt_work%bt_rem_u(i, j) = 1.0_wp
         end if
      end do
      do concurrent(j=1:ny)
         bt_work%bt_rem_u(1, j) = 1.0_wp
         bt_work%bt_rem_u(nx + 1, j) = 1.0_wp
      end do

      do concurrent(j=2:ny, i=1:nx) local(k, htot_face)
         htot_face = 0.0_wp
         do k = 1, nz
            htot_face = htot_face + 0.5_wp*(ms%h_layer(i, j - 1, k) + ms%h_layer(i, j, k))
         end do
         if (htot_face > 0.0_wp) then
            bt_work%bt_rem_v(i, j) = htot_face/(htot_face + drag_dt)
         else
            bt_work%bt_rem_v(i, j) = 1.0_wp
         end if
      end do
      do concurrent(i=1:nx)
         bt_work%bt_rem_v(i, 1) = 1.0_wp
         bt_work%bt_rem_v(i, ny + 1) = 1.0_wp
      end do
   end subroutine compute_bt_rem