set_cor_ref_velocity Subroutine

   pure subroutine set_cor_ref_velocity(grid, bt_work, ms, from_u_av, metrics, n_inner)
      !! Fill `bt_work%cor_ref_u/v` — the barotropic velocity at which
      !! `subtract_fast_cor_ref` evaluates the Coriolis/advection
      !! reference it removes from the substep forcing (MOM6
      !! `ubt_Cor`/`vbt_Cor`).
      !!
      !! The reference MUST be the depth mean of the same layer
      !! velocity whose Coriolis-advection tendency (`cor%pv_flux_*`)
      !! was depth-averaged into `F_bt_u/v`, under the same weights.
      !! Otherwise the two do not cancel at τ=0 and the residual
      !! `f × (v̄_ref − v̄_slow)` enters EVERY barotropic substep as a
      !! near-constant forcing.  In a closed rotating basin that
      !! residual projects onto the gravest Poincaré seiche and pumps
      !! it exponentially (e-folding ~0.6 d on a 240 km f-plane square
      !! at dt = 300 s; growth rate ∝ dt and rising with `n_inner` —
      !! the fingerprint of a fixed per-substep forcing, not an inner
      !! loop instability).
      !!
      !! * `from_u_av = .false.` (`ssp_rk2`) — the slow tendencies were
      !!   evaluated on the prognostic `u^n`, whose depth mean is the
      !!   stage-entry `bt_ubt/bt_vbt` from `derive_bt_from_layers`.
      !!   A plain copy, so the arithmetic downstream is bit-identical
      !!   to reading `bt_ubt/bt_vbt` directly.
      !! * `from_u_av = .true.` (`pred_corr`) — the slow tendencies were
      !!   evaluated on the time-mean `u_av/v_av` (`run_stage_split`
      !!   step 2), so take ITS depth mean, weighted exactly as the
      !!   forcing depth-mean was (h, or h·visc_rem when `&ocean_bt_nml
      !!   forcing_visc_rem`).  MOM6 does the same by construction:
      !!   `ubt_Cor = Σ_k wt_u·U_Cor` with `U_Cor = u_av`, the velocity
      !!   its `CorAdCalc` used.
      !!
      !! **Under `&vcoord_nml zfixed_closed_faces` "the same weights" is
      !! load-bearing and was, for one release, wrong here.**  A closed
      !! layer carries exactly zero velocity (`mask_layer_velocities`)
      !! but a non-zero `h_face`, so with
      !! `φ = Σ_k h_face·open / Σ_k h_face` a FULL-column mean of `u_av`
      !! returns `φ·ū_open`, not `ū_open`.  The fast loop meanwhile
      !! integrates its live `(ζ+f)·v̄ − ∇KE` on `bt_ubt = ū_open`, so
      !! `subtract_fast_cor_ref` would remove `f·φ·v̄` where it must
      !! remove `f·v̄`, leaving `Δa_u = +f·(1−φ)·v̄` forcing EVERY
      !! barotropic substep *proportionally to the barotropic velocity* —
      !! an amplifier, not a seed.  Measured on ISOMIP+ Ocean0 (melt off,
      !! 30 d): `En` `1.295E-06 → 5.711E-08` and barotropic KE ×195
      !! smaller once `metrics` is passed, matching the `ssp_rk2` twin
      !! (whose `.false.` branch below is a plain copy of `bt_ubt`, and so
      !! could never have the defect) to 2.7 %.  `φ` is O(0.5) on a
      !! partial-step face, not O(1 − 1e-4).
      type(hgrid_t), intent(in) :: grid
      type(barotropic_workstate_t), intent(inout) :: bt_work
      type(multilayer_state_t), intent(in) :: ms
      logical, intent(in) :: from_u_av
         !! `.true.` under `split_scheme = "pred_corr"`.
      type(ocean_metrics_t), intent(in) :: metrics
         !! The closed-face mask carrier, REQUIRED — see the paragraph
         !! above, and `face_depth_mean_u`'s own `metrics` docstring.
         !! Knob off ⇒ the depth means take their original branch and
         !! this is byte-identical.
      integer, intent(in), optional :: n_inner
         !! Barotropic substep count, forwarded to `face_depth_mean_rem_u/v`'s
         !! MOM6 `wt_u` floor when `bt_forcing_visc_rem` is on.  Optional
         !! (defaults to 1) ONLY so call sites that never set
         !! `forcing_visc_rem` (that branch is then never taken) need not
         !! be touched — a real `forcing_visc_rem` run must pass the true
         !! value or the floor's `Instep` is wrong.

      integer :: i, j, nu, nv, nx, ny, n_inner_use
      logical :: use_av

      n_inner_use = 1
      if (present(n_inner)) n_inner_use = n_inner
      use_av = from_u_av
      if (use_av) use_av = allocated(ms%u_av_layer) .and. allocated(ms%v_av_layer)

      if (use_av) then
         if (bt_work%bt_forcing_visc_rem) then
            call face_depth_mean_rem_u(grid, ms%u_av_layer, ms%h_layer, &
                                       bt_work%visc_rem_u, bt_work%cor_ref_u, ms%nz_ml, metrics, &
                                       n_inner_use, bt_work%bt_H_ref, bt_work%frhat_scheme)
            call face_depth_mean_rem_v(grid, ms%v_av_layer, ms%h_layer, &
                                       bt_work%visc_rem_v, bt_work%cor_ref_v, ms%nz_ml, metrics, &
                                       n_inner_use, bt_work%bt_H_ref, bt_work%frhat_scheme)
         else
            call face_depth_mean_u(grid, ms%u_av_layer, ms%h_layer, bt_work%cor_ref_u, &
                                   ms%nz_ml, metrics, bt_work%bt_H_ref, bt_work%frhat_scheme)
            call face_depth_mean_v(grid, ms%v_av_layer, ms%h_layer, bt_work%cor_ref_v, &
                                   ms%nz_ml, metrics, bt_work%bt_H_ref, bt_work%frhat_scheme)
         end if
      else
         nu = size(bt_work%bt_ubt, 1)
         ny = size(bt_work%bt_ubt, 2)
         nx = size(bt_work%bt_vbt, 1)
         nv = size(bt_work%bt_vbt, 2)
         do concurrent(j=1:ny, i=1:nu)
            bt_work%cor_ref_u(i, j) = bt_work%bt_ubt(i, j)
         end do
         do concurrent(j=1:nv, i=1:nx)
            bt_work%cor_ref_v(i, j) = bt_work%bt_vbt(i, j)
         end do
      end if
   end subroutine set_cor_ref_velocity