renormalise_meridional_flux_to_vhbt Subroutine

   pure subroutine renormalise_meridional_flux_to_vhbt(grid, metrics, this, ms, vhbt, dt, skip_walls, &
                                                       has_south, has_north, visc_rem, v_cor, &
                                                       use_por, por, use_open, open_f)
      !! Apply a uniform per-face velocity correction so
      !! `Σ_k mass_flux_y_layer(i, j, k) = vhbt(i, j)` at every face.
      !! Mirror of `renormalise_zonal_flux_to_uhbt`.  See that routine
      !! for the `skip_walls` and `has_*` (MPI seam, O4 fix) semantics.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(continuity_t), intent(in) :: this
         !! Read-only here (only the PPM edge buffers are consulted).
         !! `intent(in)` is load-bearing, not tidiness: the knob-off call
         !! sites pass one of THESE buffers as the inert `por` stand-in,
         !! and `intent(in)` turns "the callee never defines it" from a
         !! comment into a compiler-enforced invariant, so the argument
         !! association can never become aliasing.
      type(multilayer_state_t), intent(inout) :: ms
      ! assumed-shape-ok: vhbt is a face-sized array (nx, ny+1); explicit-shape
      ! would require a separate nv=ny+1 argument that the caller doesn't pass.
      real(wp), intent(in) :: vhbt(:, :)
      real(wp), intent(in) :: dt
         !! Outer-step dt (s), for the CFL bracket on `dv`.
      logical, intent(in), optional :: skip_walls
      logical, intent(in), optional :: has_south, has_north
         !! Physical-edge flags (default .true. = single-rank behaviour,
         !! bit-identical); .false. at a y-decomposition seam forces the
         !! local wall-position face to be renormalised as interior.
      ! assumed-shape-ok: mirrors the `vhbt` waiver above.
      real(wp), intent(in), optional :: visc_rem(:, :, :)
         !! Per-layer viscous remnant γ_k. Absent ⇒ γ ≡ 1, bit-identical.
         !! See `renormalise_zonal_flux_to_uhbt` for the full rationale.
      ! assumed-shape-ok: face array, same waiver as `vhbt`.
      real(wp), intent(inout), optional :: v_cor(:, :, :)
         !! MOM6 `v_cor`.  Separate time-mean field, NEVER the prognostic —
         !! see the zonal twin and docs/MOM6_SPLIT_RK2_SPEC.md §5 trap 1.

      logical, intent(in) :: use_por
         !! Porous barriers active.  `.false.` => `por` is never indexed and
         !! the arithmetic below stays byte-identical to the un-narrowed form.
      real(wp), intent(in) :: por(grid%nx_total, grid%ny_total + 1, ms%nz_ml)
         !! Layer-averaged open-area fraction at this stagger (nondim),
         !! read ONLY when `use_por`.
         !!
         !! `use_por = .false.` callers must still pass a face-sized array
         !! that is genuinely device-present, and NOT the knob-off
         !! `(1,1,1)` placeholder: nvfortran builds the `do concurrent`
         !! data clause from the LOOP BOUNDS, not from the descriptor, so a
         !! placeholder is reported "partially present" and aborts under
         !! `mem:separate` even though the branch that indexes it is never
         !! taken.  The call sites therefore hand over one of this
         !! routine's own read-only PPM edge buffers as an inert stand-in
         !! (right shape, already mapped, never DEFINED here, so no
         !! argument aliasing) -- which keeps the two full-size open-area
         !! fields off the allocation list entirely for a default run.

      logical, intent(in) :: use_open
         !! z-level closed faces active.  `.false.` => `open_f` is never
         !! indexed; byte-identical to the un-masked form.  See the zonal
         !! twin for the full rationale.
      real(wp), intent(in) :: open_f(grid%nx_total, grid%ny_total + 1, ms%nz_ml)
         !! Per-layer 0/1 face-open mask at this stagger, read ONLY when
         !! `use_open`; enters the SAME weight `wk = dx_cv * por * open`.
         !! Knob-off callers pass `h_face_right_y` as the inert stand-in
         !! (never the `(1,1,1)` placeholder).

      integer :: i, j, k, nx, ny, nz, iter
      real(wp) :: v, h_face, sum_flux, sum_h, dv, target, w, wk
      real(wp) :: flux0(NZ_STACK_MAX), v0(NZ_STACK_MAX)
      real(wp) :: v_lim, dv_hi, dv_lo, vr_k, dv_k
      real(wp) :: b_lo, b_hi, dv_new
      integer :: maxit
      logical :: consistent
      logical :: skip_w, has_s, has_n, use_vr, upd_v

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      skip_w = .true.
      if (present(skip_walls)) skip_w = skip_walls
      has_s = .true.
      if (present(has_south)) has_s = has_south
      has_n = .true.
      if (present(has_north)) has_n = has_north
      upd_v = present(v_cor)
      ! See the zonal twin: gamma-weighting and v_cor act together.
      use_vr = present(visc_rem) .and. upd_v
      consistent = this%renorm_consistent_flux
      maxit = RENORM_MAXIT
      if (consistent) maxit = RENORM_MAXIT_CONSISTENT

      do concurrent(j=2:ny, i=1:nx) &
         local(k, iter, v, h_face, sum_flux, sum_h, dv, target, w, wk, flux0, v0, v_lim, dv_hi, dv_lo, &
               vr_k, dv_k, b_lo, b_hi, dv_new)
         if (skip_w .and. &
             ((j == grid%nghost + 1 .and. has_s) .or. &
              (j == grid%nghost + grid%ny_phys + 1 .and. has_n))) cycle
         w = metrics%dx_cv(i, j)
         target = vhbt(i, j)
         do k = 1, nz
            v0(k) = ms%v_face_y_layer(i, j, k)
            flux0(k) = ms%mass_flux_y_layer(i, j, k)
         end do
         ! Legacy single-step form (wet/dry composition — mirror of the
         ! zonal branch; see `renorm_legacy_single_step`).
         if (this%renorm_legacy_single_step) then
            sum_flux = 0.0_wp
            sum_h = 0.0_wp
            do k = 1, nz
               wk = w
               if (use_por) wk = w*por(i, j, k)
               if (use_open) wk = wk*open_f(i, j, k)
               if (v0(k) >= 0.0_wp) then
                  h_face = this%h_face_left_y%data(i, j, k)
               else
                  h_face = this%h_face_right_y%data(i, j, k)
               end if
               sum_flux = sum_flux + flux0(k)
               sum_h = sum_h + h_face*wk
            end do
            if (sum_h > 0.0_wp) then
               dv = (target - sum_flux)/sum_h
               do k = 1, nz
                  wk = w
                  if (use_por) wk = w*por(i, j, k)
                  if (use_open) wk = wk*open_f(i, j, k)
                  if (v0(k) >= 0.0_wp) then
                     h_face = this%h_face_left_y%data(i, j, k)
                  else
                     h_face = this%h_face_right_y%data(i, j, k)
                  end if
                  ms%mass_flux_y_layer(i, j, k) = flux0(k) + dv*h_face*wk
               end do
            end if
            cycle
         end if
         ! CFL bracket, mirror of the zonal routine — including its
         ! land-face 1/0 guard: `idyCv` is zeroed at `wet_v == 0` by
         ! `metrics_apply_land_mask`, exactly as `idxCu` is at `wet_u == 0`.
         v_lim = RENORM_CFL/max(dt*metrics%idyCv(i, j), H_DIV_EPS)
         dv_hi = huge(1.0_wp)
         dv_lo = -huge(1.0_wp)
         do k = 1, nz
            if (use_vr) then
               vr_k = visc_rem(i, j, k)
               if (vr_k > RENORM_VR_MIN) then
                  dv_hi = min(dv_hi, (v_lim - v0(k))/vr_k)
                  dv_lo = max(dv_lo, (-v_lim - v0(k))/vr_k)
               end if
            else
               dv_hi = min(dv_hi, v_lim - v0(k))
               dv_lo = max(dv_lo, -v_lim - v0(k))
            end if
         end do
         dv_hi = max(dv_hi, 0.0_wp)
         dv_lo = min(dv_lo, 0.0_wp)
         ! Newton on `dv` with the DONOR RE-PICKED each iteration (see the
         ! RENORM_MAXIT docstring; mirror of the zonal routine).
         dv = 0.0_wp
         sum_h = 0.0_wp
         b_lo = -huge(1.0_wp)
         b_hi = huge(1.0_wp)
         do iter = 1, maxit
            sum_flux = 0.0_wp
            sum_h = 0.0_wp
            do k = 1, nz
               wk = w
               if (use_por) wk = w*por(i, j, k)
               if (use_open) wk = wk*open_f(i, j, k)
               dv_k = dv
               if (use_vr) dv_k = dv*visc_rem(i, j, k)
               if (v0(k) + dv_k >= 0.0_wp) then
                  h_face = this%h_face_left_y%data(i, j, k)
               else
                  h_face = this%h_face_right_y%data(i, j, k)
               end if
               ! `consistent`: a layer whose donor FLIPPED under the
               ! correction carries the flux of its corrected velocity
               ! through its NEW donor, `(v0+dv_k)·h_face·wk` — continuous
               ! (→ 0 from both sides) where the historical
               ! `flux0 + dv_k·h_face` jumps by `v0·(h_new − h_old)·wk`.
               ! Unflipped layers keep the historical expression.
               if (consistent .and. ((v0(k) + dv_k >= 0.0_wp) .neqv. (v0(k) >= 0.0_wp))) then
                  sum_flux = sum_flux + (v0(k) + dv_k)*h_face*wk
               else
                  sum_flux = sum_flux + (flux0(k) + dv_k*h_face*wk)
               end if
               if (use_vr) then
                  sum_h = sum_h + visc_rem(i, j, k)*h_face*wk
               else
                  sum_h = sum_h + h_face*wk
               end if
            end do
            if (sum_h <= 0.0_wp) exit
            if (abs(target - sum_flux) <= RENORM_TOL*max(1.0_wp, abs(target))) exit
            if (consistent) then
               ! Monotone, continuous F(dv): keep a bracket around the root
               ! and bisect whenever the Newton step leaves it (MOM6
               ! `zonal_flux_adjust`: Newton + bisection).
               if (sum_flux < target) then
                  b_lo = max(b_lo, dv)
               else
                  b_hi = min(b_hi, dv)
               end if
               dv_new = dv + (target - sum_flux)/sum_h
               if (dv_new <= b_lo .or. dv_new >= b_hi) then
                  if (b_lo > -huge(1.0_wp) .and. b_hi < huge(1.0_wp)) dv_new = 0.5_wp*(b_lo + b_hi)
               end if
               dv = min(max(dv_new, dv_lo), dv_hi)
            else
               dv = min(max(dv + (target - sum_flux)/sum_h, dv_lo), dv_hi)
            end if
         end do
         if (sum_h > 0.0_wp) then
            do k = 1, nz
               wk = w
               if (use_por) wk = w*por(i, j, k)
               if (use_open) wk = wk*open_f(i, j, k)
               dv_k = dv
               if (use_vr) dv_k = dv*visc_rem(i, j, k)
               if (v0(k) + dv_k >= 0.0_wp) then
                  h_face = this%h_face_left_y%data(i, j, k)
               else
                  h_face = this%h_face_right_y%data(i, j, k)
               end if
               if (consistent .and. ((v0(k) + dv_k >= 0.0_wp) .neqv. (v0(k) >= 0.0_wp))) then
                  ms%mass_flux_y_layer(i, j, k) = (v0(k) + dv_k)*h_face*wk
               else
                  ms%mass_flux_y_layer(i, j, k) = flux0(k) + dv_k*h_face*wk
               end if
               ! MOM6 `v_cor = v + dv·visc_rem` — see the zonal twin.
               if (upd_v) v_cor(i, j, k) = v0(k) + dv_k
            end do
         end if
      end do
   end subroutine renormalise_meridional_flux_to_vhbt