compute_bt_rem_from_visc_rem Subroutine

   pure subroutine compute_bt_rem_from_visc_rem(grid, bt_work, ms, metrics, n_inner)
      !! PR-2 (bt-rem-from-av-rem): build `bt_rem_u/v` from the SAME
      !! viscous remnant the layered momentum solve uses, MOM6's
      !! barotropic solver.  Dispatched the same way as
      !! `compute_bt_rem` — a RESETTER, mutually exclusive at configure
      !! with `bt_substep_drag` (D2, double-counted bed drag) and with
      !! `bt_halo > 0` (`validate_config`) — so this and `compute_bt_rem`/
      !! `reset_bt_rem` never both run for the same stage; `src/core/
      !! ocean/README.md`'s "exactly one resets, everything else
      !! MULTIPLIES" contract gets this as its third resetter.
      !!
      !! Two steps:
      !!
      !! 1. `av_rem_u/v := Σ_k frhat_k·visc_rem_k`, `frhat_k` the PLAIN
      !!    face-thickness fraction (`h_face_k / Σ_k h_face_k`) —
      !!    **not** `face_depth_mean_rem_u`'s `wt_u = h_face·visc_rem`
      !!    weighting (that one is MOM6's FORCING weight, `forcing_
      !!    visc_rem`/PR-3 scope; this is the plain depth mean MOM6 calls
      !!    `frhatu`). `frhat_k/Σ_k h_face_k` is EXACTLY
      !!    `face_depth_mean_u`'s own weight (num = Σ F·h_face, denom =
      !!    Σ h_face), so `av_rem_u = face_depth_mean_u(visc_rem_u,
      !!    h_layer)` — no separate kernel needed; this reuses the SAME
      !!    `h_face`/`metrics%open_u` branches `derive_bt_from_layers`
      !!    builds `ubt` with (the `metrics` REQUIRED-argument contract:
      !!    see that routine's docstring), so `av_rem` is the depth mean
      !!    over the SAME column the fast loop actually transports on.
      !!    `face_depth_mean_u` already returns `0` on a dry/fully-closed
      !!    face (`denom <= 0`), which is exactly the MOM6 "av_rem = 0 on
      !!    a massless column" edge case.
      !!
      !!    NOTE this `frhat` is roundabout's own: the two-abutting-cell
      !!    arithmetic-mean `h_face` `face_depth_mean_u`/`derive_bt_from_
      !!    layers`/`apply_bt_correction` already share, which is what
      !!    SELF-CONSISTENCY across the BT chain requires here — not
      !!    necessarily MOM6's own `frhatu`, which comes from `BT_cont`'s
      !!    face thicknesses (a different, flux-bounded construction).
      !!    Auditing that parity (or documenting the deliberate
      !!    divergence) is PR-3 scope, not this one.
      !!
      !! 2. `bt_rem = av_rem**(1/n_inner)` where `av_rem > 0` (MOM6
      !!    `Instep = 1/nstep`), else `0` — no `max(..., eps)` floor
      !!    substitute (CLAUDE.md: the thin-cell floor is the `av_rem >
      !!    0` MASK itself, ported exactly).  `bt_strong_drag` (MOM6
      !!    `BT_STRONG_DRAG`) swaps in the rational approximation
      !!    `n_inner·av_rem/(1+(n_inner-1)·av_rem)` instead.  Land/dry
      !!    faces are left to the existing `mask_bt_rem` call that always
      !!    runs last in the dispatch (same posture as `compute_bt_rem`,
      !!    which also does not self-mask) — MOM6's own `mask2dCu`
      !!    multiply is therefore redundant with, not additional to, that
      !!    final mask pass.
      !!
      !! Ghosts: both steps run over the FULL face extent (`size(...,1)`)
      !! including ghost columns/rows, matching `face_depth_mean_u`'s own
      !! convention — `visc_rem_u/v`'s ghosts are halo-valid after PR-1's
      !! `visc_rem_halo_refresh`, so `av_rem`/`bt_rem` are correct on
      !! every face the substep loop reads, not just the owned interior.
      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`/`face_depth_mean_u`.
      integer, intent(in) :: n_inner
         !! Barotropic substeps per outer step (MOM6 `nstep`; must be
         !! >= 1 — `auto_n_inner`/the namelist floor already enforce
         !! that).  `Instep = 1/n_inner`.

      integer :: i, j, nu, nv, ny_u, nx_v
      real(wp) :: instep
      real(wp) :: rn
      integer :: av_rem_scheme
         !! `frhat_scheme` is gated OFF (forced to FRHAT_ARITHMETIC) here
         !! specifically under OPEN z_fixed/zstar steps
         !! (`.not. metrics%use_closed_faces`): measured on the 1-degree
         !! Southern Ocean open-step case, HYBRID weighting of `av_rem`
         !! suppresses the glue-damped THIN layer's (low `visc_rem`)
         !! contribution MORE than its arithmetic share, which pulls
         !! `av_rem` UP (less barotropic damping, the opposite of the
         !! intended day-253/305 fix) — hand-verified on the 188 m/8.4 m
         !! sill fixture: av_rem = 0.495 (hybrid) vs 0.379 (arithmetic) at
         !! the sill face. Under CLOSED faces this is masked downstream
         !! (`metrics%open_u/open_v` zero the filler-adjacent weight
         !! regardless of its value) and HYBRID `av_rem` is required for
         !! the closed-faces day-305 fix (`docs/visc_rem_bt_rem_plan.md`
         !! Section 7); under OPEN faces there is no such mask and the
         !! same sign flip compounds into runaway barotropic growth (En
         !! 8x the arithmetic baseline by day 9 of a 1-degree SO open-step
         !! run, still climbing) -- confirmed by bisection against every
         !! OTHER frhat call site (`derive_bt_from_layers`,
         !! `face_depth_mean_u/v`'s slow/fast forcing, `set_cor_ref_
         !! velocity`, `apply_bt_correction`'s folds), all of which stay
         !! healthy under `frhat_scheme = "hybrid"` on the SAME open-step
         !! case. Root mechanism not fully closed out (why the reference
         !! diagnostic's harmonic av_rem reportedly helped day-305 while
         !! this port's MOM6-faithful HYBRID sweep has the opposite sign
         !! is an open question for the maintainer); this gate is the
         !! conservative fix that keeps both acceptance cases healthy
         !! without re-deriving that mechanism under time pressure.

      av_rem_scheme = merge(bt_work%frhat_scheme, FRHAT_ARITHMETIC, metrics%use_closed_faces)
      call face_depth_mean_u(grid, bt_work%visc_rem_u, ms%h_layer, bt_work%av_rem_u, &
                             ms%nz_ml, metrics, bt_work%bt_H_ref, av_rem_scheme)
      call face_depth_mean_v(grid, bt_work%visc_rem_v, ms%h_layer, bt_work%av_rem_v, &
                             ms%nz_ml, metrics, bt_work%bt_H_ref, av_rem_scheme)

      nu = size(bt_work%av_rem_u, 1)
      ny_u = size(bt_work%av_rem_u, 2)
      nx_v = size(bt_work%av_rem_v, 1)
      nv = size(bt_work%av_rem_v, 2)
      instep = 1.0_wp/real(n_inner, wp)
      rn = real(n_inner, wp)

      if (bt_work%bt_strong_drag) then
         do concurrent(j=1:ny_u, i=1:nu)
            bt_work%bt_rem_u(i, j) = 0.0_wp
            if (bt_work%av_rem_u(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_u(i, j))) then
               bt_work%bt_rem_u(i, j) = (rn*bt_work%av_rem_u(i, j))/ &
                                        (1.0_wp + (rn - 1.0_wp)*bt_work%av_rem_u(i, j))
            end if
         end do
         do concurrent(j=1:nv, i=1:nx_v)
            bt_work%bt_rem_v(i, j) = 0.0_wp
            if (bt_work%av_rem_v(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_v(i, j))) then
               bt_work%bt_rem_v(i, j) = (rn*bt_work%av_rem_v(i, j))/ &
                                        (1.0_wp + (rn - 1.0_wp)*bt_work%av_rem_v(i, j))
            end if
         end do
      else
         do concurrent(j=1:ny_u, i=1:nu)
            bt_work%bt_rem_u(i, j) = 0.0_wp
            if (bt_work%av_rem_u(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_u(i, j))) then
               bt_work%bt_rem_u(i, j) = bt_work%av_rem_u(i, j)**instep
            end if
         end do
         do concurrent(j=1:nv, i=1:nx_v)
            bt_work%bt_rem_v(i, j) = 0.0_wp
            if (bt_work%av_rem_v(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_v(i, j))) then
               bt_work%bt_rem_v(i, j) = bt_work%av_rem_v(i, j)**instep
            end if
         end do
      end if
   end subroutine compute_bt_rem_from_visc_rem