apply_bt_correction Subroutine

   pure subroutine apply_bt_correction(bt_work, ms, dt, metrics, skip_h_rescale, &
                                       grid, use_bc_pgf, use_visc_rem, scale, n_nonfin, n_inner)
      !! Replace the bt mode in the per-layer face velocities with the
      !! barotropic-substep end-step value, adding `Δu·wt_k` to every layer,
      !! `Δu = u_bt_end − u_bt_at_n − dt·F_bt_u` (same for v). Split-explicit
      !! convention (Hallberg 2009): momentum uses the END-of-step barotropic
      !! velocity; layer continuity earlier used the time-mean transports.
      !! Both legs of the corrector use the same end-step anchor (mismatched
      !! anchors overshoot the gravity-wave phase speed). Also rescales
      !! `h_layer` uniformly so the column total matches `H_ref + η_end`.
      !! `hTr` is deliberately NOT rescaled (would break exact tracer mass
      !! conservation; T = hTr/h drifts by O((η_end−η*_slow)/H) per step).
      !!
      !! ### The fold weight
      !!
      !! `wt_k = open_k·vr_k / ⟨vr⟩_h`,  `⟨vr⟩_h = Σ_k h_o·vr_k / Σ_k h_o`,
      !! `h_o = h_face·open`, so `Σ_k h_o·wt_k = Σ_k h_o` and the OPEN-column
      !! depth mean moves by exactly `Δu`.  `vr_k = visc_rem(k)` with
      !! `use_visc_rem` (MOM6 `visc_rem_u`; the drag-aware weighting),
      !! else 1, and `open ≡ 1` unless `metrics%use_closed_faces`.  With
      !! neither it is the uniform fold, `wt ≡ 1` — MOM6's own barotropic
      !! acceleration, `accel_layer_u(I,j,k) = u_accel_bt(I,j)`, the
      !! same increment in every layer.
      !!
      !! The weight does NOT carry `h`.  An earlier opt-in h-weighted form
      !! (`wt = h_face/⟨h⟩_h`, `&ocean_bt_nml correction_h_weighted`, now
      !! retired and refused by `validate_config`) adds, beyond the
      !! barotropic `ΔKE`, a positive-definite source `½Δ²·H·(κ−1)`,
      !! `κ = Σh³Σh/(Σh²)² ≥ 1`, plus a shear feedback `Δ·Σ h·u′·(wt−1)`
      !! that the fold keeps feeding: on a stretched `z_fixed` stack
      !! (`κ−1 ≈ 0.2`) the feedback ran 15-28x the source and grew the
      !! 1-degree Southern Ocean and the coastal-noise box ~5-8x until
      !! they went non-finite.  `frhatu·visc_rem` is MOM6's AVERAGING
      !! weight (BT_force, ubt), never its distribution.
      !!
      !! `skip_h_rescale` — disable the h-rescale (Lagrangian vcoord, where slow
      !!   continuity's Σh_layer is authoritative and the ALE remap relayers).
      !!   Default `.false.`.
      !! `use_bc_pgf` — add the per-layer baroclinic-PGF retro-correction
      !!   Δu_bc = -dt·((pbce(R,k)-gtot_W(R))·e_anom(R) -
      !!   (pbce(L,k)-gtot_E(L))·e_anom(L))/dx. Depth-mean zero by construction,
      !!   so the mass-flux invariant survives. Requires `grid` present
      !!   and pbce/gtot_*/e_anom populated by the caller. No-op when omitted.
      type(barotropic_workstate_t), intent(in) :: bt_work
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      type(ocean_metrics_t), intent(in) :: metrics
         !! Curvilinear horizontal metrics — read by the bc-PGF
         !! retro-correction (`use_bc_pgf = .true.` divides the e_anom
         !! gradient by `idxCu`/`idyCv`), and the carrier of the z-level
         !! closed-face mask.
         !!
         !! **REQUIRED**, like the same argument of
         !! `derive_bt_from_layers`, `face_depth_mean_u/v` and
         !! `set_cor_ref_velocity`: an OPTIONAL `metrics` that silently
         !! selects the full-column fold when omitted is the defect class
         !! that shipped the `set_cor_ref_velocity` bug.  A caller that
         !! wants the full-column fold passes a metrics object whose
         !! `use_closed_faces` latch is `.false.` (the default).
         !!
         !! When `metrics%use_closed_faces` is set a CLOSED layer gets
         !! `wt = 0` and receives nothing, BY CONSTRUCTION — the weight
         !! carries `open` and `⟨vr⟩_h` is the OPEN-column mean, so the
         !! OPEN-column depth mean (the column `derive_bt_from_layers` and
         !! `face_depth_mean_*` weight by once the knob is on) shifts by
         !! exactly `Δu`.  This replaces the spike's fold-then-mask: there
         !! `Δu` was added uniformly to every layer and
         !! `mask_layer_velocities` removed it again from the closed ones,
         !! so the layer depth mean fell short of `ubt_end` by
         !! `Δu·(Σ_closed h)/(Σ_k h)` — preserved by CANCELLATION rather
         !! than by construction.
      logical, intent(in), optional :: skip_h_rescale
      type(hgrid_t), intent(in), optional :: grid
      logical, intent(in), optional :: use_bc_pgf
      logical, intent(in), optional :: use_visc_rem
         !! Weight the fold by `visc_rem/⟨visc_rem⟩_h`.  RETIRED as a live
         !! namelist path (`&ocean_bt_nml correction_visc_rem` is
         !! fail-loud at configure, D1 follow-up — MOM6's `accel_layer_u`
         !! never weights this fold) but the kernel dispatch stays, for
         !! its own direct unit tests.  Default `.false.` ⇒ the uniform
         !! fold, which `visc_rem_chain` also uses.
      real(wp), intent(in), optional :: scale
         !! Multiplier on the Δu correction (default 1, bit-identical).
         !! The pred_corr PREDICTOR passes `BE` so the provisional velocity
         !! is `up = u + dt_pred·(u_bc_accel + u_accel_bt)` with
         !! `dt_pred = BE·dt` (SPEC §2 P8) — the tendency applies are
         !! scaled by BE at their call sites, and this scales the
         !! barotropic-increment leg to match.
      integer, intent(out), optional :: n_nonfin
         !! Count of faces whose barotropic-correction Δ (`bt_*_end − *_at_n −
         !! dt·F_bt`) is NON-FINITE.  In a supercritical hot state the BT
         !! substep loop can reach Inf on at-floor columns, and the fold's
         !! `finite − Inf` mints NaN into the layer velocity; the guard SKIPS
         !! the fold write for such a face (leaving its velocity as-is for the
         !! truncation's NaN-catch backstop) and this counts it loudly.  0 on
         !! a healthy run.
      integer, intent(in), optional :: n_inner
         !! Barotropic substeps per outer step.  REQUIRED when
         !! `bt_work%bt_rescale_strong_drag` is on (PR-2, MOM6
         !! `RESCALE_STRONG_DRAG`):
         !! `bt_strong_drag`'s rational-approximation `bt_rem` does not
         !! satisfy `bt_rem**n_in == av_rem` exactly (unlike the plain
         !! power form, which does by construction), so the Δu/Δv
         !! correction is rescaled by `min(bt_rem**n_in/av_rem, 1.0)`
         !! before being distributed into the layers — keeping the
         !! correction consistent with the TRUE depth-mean remnant.
         !! Ignored when `bt_rescale_strong_drag` is off.

      integer :: i, j, k, nu, nv, nx, ny, nz, nfin
      integer :: n_in
         !! Local copy of the optional `n_inner`.  Never read an absent
         !! optional inside `do concurrent`: ifx lowers the loop to an
         !! OpenMP region and loads every captured scalar at region entry,
         !! so a null reference segfaults even on a branch never taken.
      real(wp) :: delta_u, delta_v, total_h_old, total_h_new, ratio
      real(wp) :: du_scale
      real(wp) :: h_face, sum_h, sum_hvr, vr_bar, wt, vr_k
      real(wp) :: du_bc, dv_bc
      real(wp) :: href_l, href_r, e_prev
      real(wp) :: hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc
         !! Hand-inlined `frhat_h_face_step` locals for the visc_rem-weighted
         !! fold's two `do concurrent` blocks below -- see the comment at
         !! their first use for why this is inlined rather than called.
      integer :: il, ir, jl, jr, frhat_scheme
      logical :: do_rescale, do_bc_pgf, do_visc_rem, do_open, do_bt_rescale

      ! FRHAT_HYBRID only under closed faces -- see `derive_bt_from_layers`'
      ! matching comment (the same gate `compute_bt_rem_from_visc_rem`
      ! already applies to av_rem, extended here to this fold's two
      ! `do concurrent` blocks and their hand-inlined GPU twins below).
      frhat_scheme = merge(bt_work%frhat_scheme, FRHAT_ARITHMETIC, metrics%use_closed_faces)
      do_rescale = .true.
      if (present(skip_h_rescale)) do_rescale = .not. skip_h_rescale
      do_bc_pgf = .false.
      if (present(use_bc_pgf)) do_bc_pgf = use_bc_pgf
      do_visc_rem = .false.
      if (present(use_visc_rem)) do_visc_rem = use_visc_rem
      du_scale = 1.0_wp
      if (present(scale)) du_scale = scale
      do_open = metrics%use_closed_faces
      n_in = 1
      if (present(n_inner)) n_in = n_inner
      do_bt_rescale = bt_work%bt_rescale_strong_drag
      if (do_bt_rescale .and. .not. present(n_inner)) then
         error stop "apply_bt_correction: rescale_strong_drag requires n_inner"
      end if
      if (do_bc_pgf .and. .not. present(grid)) then
         error stop "apply_bt_correction: use_bc_pgf=.true. requires grid"
      end if

      nu = size(ms%u_face_x_layer, 1)
      ny = size(ms%u_face_x_layer, 2)
      nx = size(ms%v_face_y_layer, 1)
      nv = size(ms%v_face_y_layer, 2)
      nz = ms%nz_ml

      ! Loud count of non-finite fold Δ (2D, once per face) — the BT loop can
      ! reach Inf on at-floor columns in a supercritical state.  Read-only
      ! reduction (write out of the reduction loop, per the truncation fix).
      ! (No `present()` clause: this kernel's fold uses stdpar `do concurrent`
      ! managed memory, so callers do not acc-map bt_work; present_or_copyin
      ! reads the device copy in production and copies-in host data in the
      ! unmapped unit tests.)
      if (present(n_nonfin)) then
         nfin = 0
         do concurrent(j=1:ny, i=1:nu) reduce(+:nfin)
            if (.not. ieee_is_finite(bt_work%bt_ubt_end(i, j) - bt_work%ubt_at_n(i, j) &
                                     - dt*bt_work%F_bt_u(i, j))) nfin = nfin + 1
         end do
         do concurrent(j=1:nv, i=1:nx) reduce(+:nfin)
            if (.not. ieee_is_finite(bt_work%bt_vbt_end(i, j) - bt_work%vbt_at_n(i, j) &
                                     - dt*bt_work%F_bt_v(i, j))) nfin = nfin + 1
         end do
         n_nonfin = nfin
      end if

      if (do_open) then
         ! OPEN-LAYER fold (`&vcoord_nml zfixed_closed_faces`):
         ! `wt = open·vr/⟨vr⟩_h` over `h_o = h_face·open`.  A CLOSED layer
         ! receives nothing.  Without visc_rem `wt = open` — written as
         ! exactly that (not `open·1/1`) so the no-visc_rem closed-face
         ! answer is the one this branch always gave.
         do concurrent(j=1:ny, i=1:nu) &
            local(k, il, ir, href_l, href_r, e_prev, delta_u, sum_h, sum_hvr, h_face, vr_bar, wt, vr_k, &
                  hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
            ! `frhat_h_face_step` hand-inlined -- see the comment at the
            ! visc_rem-weighted fold below (same `apply_bt_correction`
            ! call-site miscompile under NVHPC 25.5 `-stdpar=gpu`,
            ! reproduced here too via compute-sanitizer).
            delta_u = du_scale*(bt_work%bt_ubt_end(i, j) - bt_work%ubt_at_n(i, j) - dt*bt_work%F_bt_u(i, j))
            if (do_bt_rescale) then
               if (bt_work%av_rem_u(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_u(i, j))) then
                  delta_u = delta_u*min(bt_work%bt_rem_u(i, j)**n_in/bt_work%av_rem_u(i, j), 1.0_wp)
               end if
            end if
            il = max(1, i - 1)
            ir = min(nu - 1, i)
            href_l = bt_work%bt_H_ref(il, j)
            href_r = bt_work%bt_H_ref(ir, j)
            e_prev = -0.5_wp*(href_l + href_r)
            sum_h = 0.0_wp
            sum_hvr = 0.0_wp
            do k = 1, nz
               hl_loc = ms%h_layer(il, j, k)
               hr_loc = ms%h_layer(ir, j, k)
               h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
               if (frhat_scheme == FRHAT_HYBRID) then
                  d_shallow_loc = -min(href_l, href_r)
                  e_cur_loc = e_prev + h_arith_loc
                  ! vanished-ok: hand-inlined frhat_h_face_step copy (NVHPC call-site miscompile, see local() comment above)
                  if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                     h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                  else if (e_prev >= d_shallow_loc) then
                     h_face = h_arith_loc
                  else
                     h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                     if (e_cur_loc <= d_shallow_loc) then
                        h_face = h_harm_loc
                     else
                        wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                        h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                     end if
                  end if
                  e_prev = e_cur_loc
               else
                  h_face = h_arith_loc
               end if
               h_face = h_face*metrics%open_u(i, j, k)
               if (do_visc_rem) then
                  vr_k = bt_work%visc_rem_u(i, j, k)
               else
                  vr_k = 1.0_wp
               end if
               sum_h = sum_h + h_face
               sum_hvr = sum_hvr + h_face*vr_k
            end do
            if (sum_h > 0.0_wp .and. ieee_is_finite(delta_u)) then
               vr_bar = 1.0_wp
               if (do_visc_rem .and. sum_hvr > 0.0_wp) vr_bar = sum_hvr/sum_h
               do k = 1, nz
                  if (do_visc_rem) then
                     wt = metrics%open_u(i, j, k)*bt_work%visc_rem_u(i, j, k)/vr_bar
                  else
                     wt = metrics%open_u(i, j, k)
                  end if
                  ms%u_face_x_layer(i, j, k) = ms%u_face_x_layer(i, j, k) + delta_u*wt
               end do
            end if
         end do
         do concurrent(j=1:nv, i=1:nx) &
            local(k, jl, jr, href_l, href_r, e_prev, delta_v, sum_h, sum_hvr, h_face, vr_bar, wt, vr_k, &
                  hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
            ! See the matching comment in the u-branch above.
            delta_v = du_scale*(bt_work%bt_vbt_end(i, j) - bt_work%vbt_at_n(i, j) - dt*bt_work%F_bt_v(i, j))
            if (do_bt_rescale) then
               if (bt_work%av_rem_v(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_v(i, j))) then
                  delta_v = delta_v*min(bt_work%bt_rem_v(i, j)**n_in/bt_work%av_rem_v(i, j), 1.0_wp)
               end if
            end if
            jl = max(1, j - 1)
            jr = min(nv - 1, j)
            href_l = bt_work%bt_H_ref(i, jl)
            href_r = bt_work%bt_H_ref(i, jr)
            e_prev = -0.5_wp*(href_l + href_r)
            sum_h = 0.0_wp
            sum_hvr = 0.0_wp
            do k = 1, nz
               hl_loc = ms%h_layer(i, jl, k)
               hr_loc = ms%h_layer(i, jr, k)
               h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
               if (frhat_scheme == FRHAT_HYBRID) then
                  d_shallow_loc = -min(href_l, href_r)
                  e_cur_loc = e_prev + h_arith_loc
                  ! vanished-ok: hand-inlined frhat_h_face_step copy, same reason as the u-branch above
                  if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                     h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                  else if (e_prev >= d_shallow_loc) then
                     h_face = h_arith_loc
                  else
                     h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                     if (e_cur_loc <= d_shallow_loc) then
                        h_face = h_harm_loc
                     else
                        wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                        h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                     end if
                  end if
                  e_prev = e_cur_loc
               else
                  h_face = h_arith_loc
               end if
               h_face = h_face*metrics%open_v(i, j, k)
               if (do_visc_rem) then
                  vr_k = bt_work%visc_rem_v(i, j, k)
               else
                  vr_k = 1.0_wp
               end if
               sum_h = sum_h + h_face
               sum_hvr = sum_hvr + h_face*vr_k
            end do
            if (sum_h > 0.0_wp .and. ieee_is_finite(delta_v)) then
               vr_bar = 1.0_wp
               if (do_visc_rem .and. sum_hvr > 0.0_wp) vr_bar = sum_hvr/sum_h
               do k = 1, nz
                  if (do_visc_rem) then
                     wt = metrics%open_v(i, j, k)*bt_work%visc_rem_v(i, j, k)/vr_bar
                  else
                     wt = metrics%open_v(i, j, k)
                  end if
                  ms%v_face_y_layer(i, j, k) = ms%v_face_y_layer(i, j, k) + delta_v*wt
               end do
            end if
         end do
      else if (.not. do_visc_rem) then
         ! Uniform Δu distribution — every layer gets the same Δu.  The finite
         ! guard skips a face whose Δ is non-finite (Inf/NaN from a blown-up BT
         ! loop) so the fold never mints NaN into the layer velocity.
         do concurrent(k=1:nz, j=1:ny, i=1:nu) local(delta_u)
            delta_u = du_scale*(bt_work%bt_ubt_end(i, j) - bt_work%ubt_at_n(i, j) - dt*bt_work%F_bt_u(i, j))
            if (do_bt_rescale) then
               if (bt_work%av_rem_u(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_u(i, j))) then
                  delta_u = delta_u*min(bt_work%bt_rem_u(i, j)**n_in/bt_work%av_rem_u(i, j), 1.0_wp)
               end if
            end if
            if (ieee_is_finite(delta_u)) then
               ms%u_face_x_layer(i, j, k) = ms%u_face_x_layer(i, j, k) + delta_u
            end if
         end do
         do concurrent(k=1:nz, j=1:nv, i=1:nx) local(delta_v)
            delta_v = du_scale*(bt_work%bt_vbt_end(i, j) - bt_work%vbt_at_n(i, j) - dt*bt_work%F_bt_v(i, j))
            if (do_bt_rescale) then
               if (bt_work%av_rem_v(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_v(i, j))) then
                  delta_v = delta_v*min(bt_work%bt_rem_v(i, j)**n_in/bt_work%av_rem_v(i, j), 1.0_wp)
               end if
            end if
            if (ieee_is_finite(delta_v)) then
               ms%v_face_y_layer(i, j, k) = ms%v_face_y_layer(i, j, k) + delta_v
            end if
         end do
      else
         ! visc_rem-weighted fold on the full column: `wt = vr/⟨vr⟩_h`.
         ! A column with no thickness, or with `Σ h·vr = 0` (every layer
         ! fully damped), takes the uniform increment: `wt ≡ 1` is the
         ! only weight with the right depth mean when `⟨vr⟩_h` is
         ! undefined.
         do concurrent(j=1:ny, i=1:nu) &
            local(k, il, ir, href_l, href_r, e_prev, delta_u, sum_h, sum_hvr, h_face, vr_bar, wt, &
                  hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
            ! `frhat_h_face_step` is called BY HAND here rather than via `call`
            ! (contrast `derive_bt_from_layers`/`face_depth_mean_*`, which call
            ! it directly and run clean on GPU): under NVHPC 25.5
            ! `-stdpar=gpu,mem:separate` this specific call site faulted with
            ! an "Invalid __global__ read" inside the callee (compute-
            ! sanitizer memcheck, 2026-10-06) -- a device-codegen defect, not
            ! a bounds bug (`il`/`ir` are always in `[1, nx]` by construction)
            ! and not a missing-device-mapping bug (reproduced even with
            ! `bt_work`/`ms` fully `enter_data`-mapped). Moving the `il`/`ir`
            ! assignment out of the `ieee_is_finite` guard did not fix it
            ! either. The inline copy below is textually identical to
            ! `rdb_frhat_face.inc`'s body (CLAUDE.md's "Cross-TU helper
            ! inlining" gotcha, applied one step further: inline the BODY, not
            ! just the routine, when even an in-module call to a `!$acc
            ! routine seq` helper miscompiles at a specific call site).
            il = max(1, i - 1)
            ir = min(nu - 1, i)
            href_l = bt_work%bt_H_ref(il, j)
            href_r = bt_work%bt_H_ref(ir, j)
            delta_u = du_scale*(bt_work%bt_ubt_end(i, j) - bt_work%ubt_at_n(i, j) - dt*bt_work%F_bt_u(i, j))
            if (do_bt_rescale) then
               if (bt_work%av_rem_u(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_u(i, j))) then
                  delta_u = delta_u*min(bt_work%bt_rem_u(i, j)**n_in/bt_work%av_rem_u(i, j), 1.0_wp)
               end if
            end if
            if (ieee_is_finite(delta_u)) then
               e_prev = -0.5_wp*(href_l + href_r)
               sum_h = 0.0_wp
               sum_hvr = 0.0_wp
               do k = 1, nz
                  hl_loc = ms%h_layer(il, j, k)
                  hr_loc = ms%h_layer(ir, j, k)
                  h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
                  if (frhat_scheme == FRHAT_HYBRID) then
                     d_shallow_loc = -min(href_l, href_r)
                     e_cur_loc = e_prev + h_arith_loc
                     ! vanished-ok: hand-inlined frhat_h_face_step copy (NVHPC call-site miscompile, see local() comment above)
                     if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                        h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                     else if (e_prev >= d_shallow_loc) then
                        h_face = h_arith_loc
                     else
                        h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                        if (e_cur_loc <= d_shallow_loc) then
                           h_face = h_harm_loc
                        else
                           wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                           h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                        end if
                     end if
                     e_prev = e_cur_loc
                  else
                     h_face = h_arith_loc
                  end if
                  sum_h = sum_h + h_face
                  sum_hvr = sum_hvr + h_face*bt_work%visc_rem_u(i, j, k)
               end do
               if (sum_h > 0.0_wp .and. sum_hvr > 0.0_wp) then
                  vr_bar = sum_hvr/sum_h
                  do k = 1, nz
                     wt = bt_work%visc_rem_u(i, j, k)/vr_bar
                     ms%u_face_x_layer(i, j, k) = ms%u_face_x_layer(i, j, k) + delta_u*wt
                  end do
               else
                  do k = 1, nz
                     ms%u_face_x_layer(i, j, k) = ms%u_face_x_layer(i, j, k) + delta_u
                  end do
               end if
            end if
         end do
         do concurrent(j=1:nv, i=1:nx) &
            local(k, jl, jr, href_l, href_r, e_prev, delta_v, sum_h, sum_hvr, h_face, vr_bar, wt, &
                  hl_loc, hr_loc, h_arith_loc, h_harm_loc, e_cur_loc, d_shallow_loc, wt_arith_loc)
            ! See the matching comment in the u-branch above -- `jl`/`jr`/
            ! `href_l`/`href_r` computed unconditionally, and
            ! `frhat_h_face_step` inlined by hand for the same reason.
            jl = max(1, j - 1)
            jr = min(nv - 1, j)
            href_l = bt_work%bt_H_ref(i, jl)
            href_r = bt_work%bt_H_ref(i, jr)
            delta_v = du_scale*(bt_work%bt_vbt_end(i, j) - bt_work%vbt_at_n(i, j) - dt*bt_work%F_bt_v(i, j))
            if (do_bt_rescale) then
               if (bt_work%av_rem_v(i, j) > 0.0_wp .and. ieee_is_finite(bt_work%av_rem_v(i, j))) then
                  delta_v = delta_v*min(bt_work%bt_rem_v(i, j)**n_in/bt_work%av_rem_v(i, j), 1.0_wp)
               end if
            end if
            if (ieee_is_finite(delta_v)) then
               e_prev = -0.5_wp*(href_l + href_r)
               sum_h = 0.0_wp
               sum_hvr = 0.0_wp
               do k = 1, nz
                  hl_loc = ms%h_layer(i, jl, k)
                  hr_loc = ms%h_layer(i, jr, k)
                  h_arith_loc = 0.5_wp*(hl_loc + hr_loc)
                  if (frhat_scheme == FRHAT_HYBRID) then
                     d_shallow_loc = -min(href_l, href_r)
                     e_cur_loc = e_prev + h_arith_loc
                     ! vanished-ok: hand-inlined frhat_h_face_step copy, same reason as the u-branch above
                     if (hl_loc <= H_VANISHED .or. hr_loc <= H_VANISHED) then
                        h_face = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                     else if (e_prev >= d_shallow_loc) then
                        h_face = h_arith_loc
                     else
                        h_harm_loc = (hl_loc*hr_loc)/(h_arith_loc + H_DIV_EPS)
                        if (e_cur_loc <= d_shallow_loc) then
                           h_face = h_harm_loc
                        else
                           wt_arith_loc = (e_cur_loc - d_shallow_loc)/(h_arith_loc + H_DIV_EPS)
                           h_face = wt_arith_loc*h_arith_loc + (1.0_wp - wt_arith_loc)*h_harm_loc
                        end if
                     end if
                     e_prev = e_cur_loc
                  else
                     h_face = h_arith_loc
                  end if
                  sum_h = sum_h + h_face
                  sum_hvr = sum_hvr + h_face*bt_work%visc_rem_v(i, j, k)
               end do
               if (sum_h > 0.0_wp .and. sum_hvr > 0.0_wp) then
                  vr_bar = sum_hvr/sum_h
                  do k = 1, nz
                     wt = bt_work%visc_rem_v(i, j, k)/vr_bar
                     ms%v_face_y_layer(i, j, k) = ms%v_face_y_layer(i, j, k) + delta_v*wt
                  end do
               else
                  do k = 1, nz
                     ms%v_face_y_layer(i, j, k) = ms%v_face_y_layer(i, j, k) + delta_v
                  end do
               end if
            end if
         end do
      end if

      ! bc-PGF additive per-layer correction. L = west/south cell, R =
      ! east/north cell. Each (pbce(k) - gtot_face) term is depth-mean-zero per
      ! column, so the column-mean velocity is unchanged (mass-flux invariant
      ! preserved). Interior faces only (wall faces already BT-zeroed).
      if (do_bc_pgf) then
         do concurrent(k=1:nz, j=1:ny, i=2:nu - 1) local(du_bc)
            du_bc = -dt*((bt_work%pbce(i, j, k) - bt_work%gtot_W(i, j)) &
                         *bt_work%e_anom(i, j) &
                         - (bt_work%pbce(i - 1, j, k) - bt_work%gtot_E(i - 1, j)) &
                         *bt_work%e_anom(i - 1, j))*metrics%idxCu(i, j)
            if (ieee_is_finite(du_bc)) ms%u_face_x_layer(i, j, k) = ms%u_face_x_layer(i, j, k) + du_bc
         end do
         do concurrent(k=1:nz, j=2:nv - 1, i=1:nx) local(dv_bc)
            dv_bc = -dt*((bt_work%pbce(i, j, k) - bt_work%gtot_S(i, j)) &
                         *bt_work%e_anom(i, j) &
                         - (bt_work%pbce(i, j - 1, k) - bt_work%gtot_N(i, j - 1)) &
                         *bt_work%e_anom(i, j - 1))*metrics%idyCv(i, j)
            if (ieee_is_finite(dv_bc)) ms%v_face_y_layer(i, j, k) = ms%v_face_y_layer(i, j, k) + dv_bc
         end do
      end if

      if (do_rescale) then
         do concurrent(j=1:size(ms%h_layer, 2), i=1:size(ms%h_layer, 1)) &
            local(k, total_h_old, total_h_new, ratio)
            total_h_old = 0.0_wp
            do k = 1, nz
               total_h_old = total_h_old + ms%h_layer(i, j, k)
            end do
            total_h_new = bt_work%bt_H_ref(i, j) + bt_work%bt_eta_end(i, j)
            if (total_h_old > 0.0_wp) then
               ratio = total_h_new/total_h_old
               do k = 1, nz
                  ms%h_layer(i, j, k) = ms%h_layer(i, j, k)*ratio
               end do
            end if
         end do
      end if
   end subroutine apply_bt_correction