continuity_zonal_flux Subroutine

   pure subroutine continuity_zonal_flux(grid, metrics, this, ms, dt, uhbt, bc, &
                                         visc_rem, u_cor)
      !! Zonal (x-only) PPM reconstruction + per-face mass flux on
      !! the multilayer C-grid.  Companion to
      !! `continuity_meridional_flux` for the
      !! directionally-split (Lie) continuity step.  Writes
      !! `ms%mass_flux_x_layer` and leaves `mass_flux_y_layer` /
      !! `flux_h_layer` untouched.  Wall faces at i=1 and i=nx+1
      !! are zeroed (closed-wall BC).
      !!
      !! Stencil + boundary treatment identical to the X half of
      !! `continuity_compute_fluxes` — the body was
      !! lifted verbatim and the Y block dropped.
      !!
      !! Optional `uhbt(i, j)` is the time-mean barotropic-substep transport
      !! at the east face.  When present, the per-layer mass fluxes
      !! are renormalised by a uniform velocity correction
      !! `du(i, j) = (uhbt - Σ_k mass_flux) / Σ_k h_face` so the
      !! vertical sum matches `uhbt` exactly.  This is the MOM6
      !! continuity-with-uhbt pattern.
      !! With this constraint, the slow continuity advances h_layer
      !! consistently with the barotropic-substep's η end-state — no
      !! post-rescale needed.  Absent ⇒ unconstrained PPM transport.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(continuity_t), intent(inout) :: this
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
         !! Time increment (s).  Used only for the swept-volume CFL when
         !! `this%vol_cfl = .true.`; ignored (bit-identical) otherwise.
      real(wp), intent(in), optional :: uhbt(:, :)
         !! Time-mean east-face transport from barotropic substep (m³/s),
         !! shape `(nx+1, ny)`.  Width-weighted (carries `dy_cu`) so it
         !! constrains the same transport `mass_flux_x_layer` now holds.
      type(ocean_bc_state_t), intent(in), optional :: bc
         !! Per-edge OBC tags.  Default (absent) -> OBC_WALL on both
         !! ends.  Non-WALL tags skip the wall-zero step at that edge.
      ! assumed-shape-ok: pure passthrough to the renormaliser, which carries
      ! the same waiver; never indexed here.
      real(wp), intent(in), optional :: visc_rem(:, :, :)
         !! Per-layer viscous remnant γ_k, forwarded to
         !! `renormalise_zonal_flux_to_uhbt`.  Absent ⇒ γ ≡ 1, bit-identical.
      ! assumed-shape-ok: pure passthrough to the renormaliser.
      real(wp), intent(inout), optional :: u_cor(:, :, :)
         !! Forwarded MOM6 `u_cor` destination — a SEPARATE time-mean field,
         !! never the prognostic.  Absent ⇒ flux-only, bit-identical.

      integer :: i, j, k, nx, ny, nz
      integer :: bc_w_tag, bc_e_tag
      real(wp) :: dh_m1, dh_0, dh_p1, h_left, h_right, u, h_face
      real(wp) :: hm2, hm1, h0, hp1, hp2
      real(wp) :: cfl_d, curv3_d, dh_d
      logical :: do_pos, do_volcfl, do_pd
      logical :: has_w_flux, has_e_flux, renorm_skip_walls
      real(wp) :: h_min_pos, two_h_lim

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      do_pos = this%use_ppm_limit_pos
      do_volcfl = this%vol_cfl
      h_min_pos = this%h_min
      ! P1 positive-definite reconstruction floor (MOM6's positive-definite
      ! PPM).  Scalars copied to locals so the DC loop never walks the
      ! continuity_t descriptor.  Off ⇒ untaken branch ⇒ bit-identical.
      do_pd = this%positive_definite
      two_h_lim = 2.0_wp*this%h_lim

      do concurrent(k=1:nz, j=1:ny, i=3:nx - 2) &
         local(dh_m1, dh_0, dh_p1, h_left, h_right, hm2, hm1, h0, hp1, hp2)
         h0 = ms%h_layer(i, j, k)
         hm1 = ppm_mirror_h(ms%h_layer(i - 1, j, k), h0, metrics%wet_T(i - 1, j))
         hp1 = ppm_mirror_h(ms%h_layer(i + 1, j, k), h0, metrics%wet_T(i + 1, j))
         hm2 = ppm_mirror_h(ms%h_layer(i - 2, j, k), hm1, metrics%wet_T(i - 2, j))
         hp2 = ppm_mirror_h(ms%h_layer(i + 2, j, k), hp1, metrics%wet_T(i + 2, j))
         call ppm_limited_slope(hm2, hm1, h0, dh_m1)
         call ppm_limited_slope(hm1, h0, hp1, dh_0)
         call ppm_limited_slope(h0, hp1, hp2, dh_p1)
         dh_0 = dh_0*metrics%wet_T(i - 1, j)*metrics%wet_T(i, j)*metrics%wet_T(i + 1, j)
         h_left = 0.5_wp*(hm1 + h0) - (dh_0 - dh_m1)/6.0_wp
         h_right = 0.5_wp*(h0 + hp1) - (dh_p1 - dh_0)/6.0_wp
         call ppm_cell_limiter(h0, h_left, h_right)
         if (do_pos) call ppm_limit_pos(h0, h_left, h_right, h_min_pos)
         if (do_pd) then
            h_left = max(h_left, two_h_lim)
            h_right = max(h_right, two_h_lim)
         end if
         this%h_face_right_x%data(i, j, k) = h_left
         this%h_face_left_x%data(i + 1, j, k) = h_right
      end do
      do concurrent(k=1:nz, j=1:ny)
         this%h_face_left_x%data(1, j, k) = ms%h_layer(1, j, k)
         this%h_face_right_x%data(1, j, k) = ms%h_layer(1, j, k)
         this%h_face_left_x%data(2, j, k) = ms%h_layer(1, j, k)
         this%h_face_right_x%data(2, j, k) = ms%h_layer(2, j, k)
         ! Face 3: cell 2's right edge falls back to 1st order (its
         ! 5-point stencil needs cell 0); cell 3's left edge is set
         ! by the interior loop above.
         this%h_face_left_x%data(3, j, k) = ms%h_layer(2, j, k)
         this%h_face_right_x%data(3, j, k) = ms%h_layer(2, j, k)
         ! Face nx-1: mirror of face 3.  Cell nx-1's left edge falls
         ! back to 1st order; cell nx-2's right edge came from the
         ! interior loop.
         this%h_face_right_x%data(nx - 1, j, k) = ms%h_layer(nx - 1, j, k)
         this%h_face_left_x%data(nx, j, k) = ms%h_layer(nx - 1, j, k)
         this%h_face_right_x%data(nx, j, k) = ms%h_layer(nx, j, k)
         this%h_face_left_x%data(nx + 1, j, k) = ms%h_layer(nx, j, k)
         this%h_face_right_x%data(nx + 1, j, k) = ms%h_layer(nx, j, k)
      end do

      do concurrent(k=1:nz, j=1:ny, i=2:nx) &
         local(u, h_face, cfl_d, curv3_d, dh_d)
         u = ms%u_face_x_layer(i, j, k)
         if (u >= 0.0_wp) then
            ! Donor = cell i-1; downwind (east) edge = h_face_left_x(i).
            h_face = this%h_face_left_x%data(i, j, k)
            if (do_volcfl) then
               ! Swept-oriented donor edge diff dh = h_L - h_R (west-east).
               dh_d = this%h_face_right_x%data(i - 1, j, k) - h_face
               curv3_d = this%h_face_right_x%data(i - 1, j, k) + h_face &
                         - 2.0_wp*ms%h_layer(i - 1, j, k)
               cfl_d = u*dt*metrics%dy_cu(i, j)*metrics%iareaT(i - 1, j)
               h_face = volcfl_face(h_face, dh_d, curv3_d, cfl_d)
            end if
         else
            ! Donor = cell i; downwind (west) edge = h_face_right_x(i).
            h_face = this%h_face_right_x%data(i, j, k)
            if (do_volcfl) then
               ! Swept-oriented donor edge diff dh = h_R - h_L (east-west).
               dh_d = this%h_face_left_x%data(i + 1, j, k) - h_face
               curv3_d = h_face + this%h_face_left_x%data(i + 1, j, k) &
                         - 2.0_wp*ms%h_layer(i, j, k)
               cfl_d = (-u)*dt*metrics%dy_cu(i, j)*metrics%iareaT(i, j)
               h_face = volcfl_face(h_face, dh_d, curv3_d, cfl_d)
            end if
         end if
         ms%mass_flux_x_layer(i, j, k) = u*h_face*metrics%dy_cu(i, j)
      end do
      do concurrent(k=1:nz, j=1:ny)
         ms%mass_flux_x_layer(1, j, k) = 0.0_wp
         ms%mass_flux_x_layer(nx + 1, j, k) = 0.0_wp
      end do
      ! Zero the PHYSICAL wall faces (not just the array edges).  The
      ! ocean slow-path apply kernels (surface stress / Coriolis / PGF
      ! / drag) don't enforce u_face_x = 0 at the physical interior
      ! boundary, so even with `u^n_wall = 0` the velocity applies
      ! between stages can drive a small `u^{stage1}_wall` that the
      ! next stage's `mass_flux = u·h_face` would treat as a wall
      ! flux.  Without this zeroing, the PPM scheme transports
      ! tracer mass across the physical wall into ghost cells, and
      ! `sum(hTr)` drifts by ~1e-6 over hundreds of outer steps.
      ! Mirrors MOM6's `mask2dCu` face mask convention.  Under MPI
      ! decomposition the gate now reads bc%has_west / bc%has_east
      ! (the flag mirrors decomp%has_west / has_east set at driver
      ! setup) so a subdomain seam is never hard-zeroed here — the
      ! halo exchange corrects it.
      !
      ! OBC dispatch: non-WALL edges keep the computed mass flux so
      ! the downstream transport sees the barotropic-substep's open-boundary
      ! velocity.  Caller must supply matching `uhbt` at those faces
      ! (the driver no longer zeros it for non-WALL edges).
      bc_w_tag = OBC_WALL
      bc_e_tag = OBC_WALL
      ! Physical-edge flags: default .true. => single-rank bit-identical
      ! behaviour; .false. at an MPI seam skips the hard zero (O0, D0).
      has_w_flux = .true.
      has_e_flux = .true.
      if (present(bc)) then
         bc_w_tag = ocean_bc_outer_face_tag(bc%west%bc_type)
         bc_e_tag = ocean_bc_outer_face_tag(bc%east%bc_type)
         has_w_flux = bc%has_west
         has_e_flux = bc%has_east
      end if
      do concurrent(k=1:nz, j=1:ny)
         if (bc_w_tag == OBC_WALL .and. has_w_flux) ms%mass_flux_x_layer(grid%nghost + 1, j, k) = 0.0_wp
         if (bc_e_tag == OBC_WALL .and. has_e_flux) ms%mass_flux_x_layer(grid%nghost + grid%nx_phys + 1, j, k) = 0.0_wp
      end do

      ! ---- Porous barriers (Adcroft 2013) ----
      ! Narrow the layer transport by the OPEN-AREA fraction of the face.
      ! Host-side gate: with the knob off there is no kernel launch and
      ! the loops above are textually unchanged, so the whole path is
      ! byte-identical.  Applied AFTER the wall zeroing (0 stays 0) and
      ! BEFORE the barotropic renormalisation, which must see the narrowed
      ! transports it is constraining.
      !
      ! WRITTEN INLINE, not as a call to `porous_narrow_3d`.  Handing
      ! `ms%mass_flux_*_layer` to an external subroutine as an
      ! `intent(inout)` actual makes nvfortran treat the array as ESCAPING,
      ! which pessimises every `do concurrent` in this routine — even
      ! though the branch never runs with the knob off.  Measured on a
      ! 600x600x50 default-path (porous OFF) double-gyre, V100: the call
      ! form costs `ocean_continuity` 10.58 s vs 9.11 s inline (+4.8% on
      ! total solver time vs origin/main; the inline form is +0.5%, i.e.
      ! noise).  `rdb_coriolis_adv` keeps the shared `porous_narrow_3d`
      ! helper — measured there at no cost.
      if (metrics%use_porous) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
            ms%mass_flux_x_layer(i, j, k) = ms%mass_flux_x_layer(i, j, k)* &
                                            metrics%por_face_area_u(i, j, k)
         end do
      end if

      ! ---- z-level closed faces: see the composition rule on
      ! `ocean_metrics_t%open_v`.  Also before the renormalisation — the
      ! renormaliser must distribute `uhbt` over the OPEN layers ONLY, so
      ! it takes the SAME mask as a weight (`use_open`/`open` below).
      if (metrics%use_closed_faces) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
            ms%mass_flux_x_layer(i, j, k) = ms%mass_flux_x_layer(i, j, k)* &
                                            metrics%open_u(i, j, k)
         end do
      end if

      ! ---- MOM6-style transport constraint ----
      ! Renormalise the per-layer mass flux so its vertical sum
      ! matches the barotropic-substep's time-mean transport.  Single Newton
      ! pass: assuming the upwind donor doesn't flip under the
      ! correction (true for small `du`), `Σ_k (u_k + du) · h_face_k =
      ! uhbt` solves to `du = (uhbt - Σ_k u_k·h_face_k) / Σ_k h_face_k`.
      ! Apply `du` per layer; the sign of u_face decides the donor
      ! (left vs right h_face), unchanged from above.
      !
      ! Wall faces (i=1, i=nx+1, plus the physical-wall zeroing
      ! above) keep mass_flux at zero and are skipped via uhbt(I,j)=0
      ! at those indices (caller's barotropic substep initialises and only
      ! accumulates on owned interior faces).
      if (present(uhbt)) then
         ! Any non-WALL zonal edge (periodic OR open-class Flather/tidal/
         ! Chapman/clamped) carries genuine transport at the physical-wall
         ! face: its mass flux was NOT zeroed above, and the barotropic
         ! substep supplies a nonzero uhbt there.  Those faces must be
         ! renormalised like interior faces so the per-layer fluxes sum to
         ! the BT transport — otherwise the boundary leaks mode mismatch.
         ! For a genuine closed WALL (uhbt = 0, mass_flux = 0) the wall face
         ! is skipped (default skip_walls=.true.): nothing to renormalise.
         ! All-WALL → skip_walls stays .true. → bit-identical to pre-OBC.
         ! `visc_rem` / `u_cor` forward straight through: an ABSENT
         ! optional passed to an optional dummy stays absent, so the no-knob
         ! path reaches the renormaliser exactly as before (bit-identical).
         ! ONE call site per `por` actual, not one per BC branch.  `skip_w`
         ! and the `has_*` flags are now computed here and always passed:
         ! when `skip_walls` is .false. the callee never consults `has_*`
         ! (the test is `skip_w .and. (...)`), so folding the two BC
         ! branches into one call is exact, not merely equivalent.  Keeping
         ! them separate cost a second inlined copy of the renormaliser's
         ! `do concurrent` per branch, which is where the measured
         ! default-path regression lived.
         !
         ! The `por` actual still differs per branch because a knob-off run
         ! must NOT hand over the `(1,1,1)` placeholder — see the `por`
         ! dummy's docstring.  `h_face_left_x` is the inert stand-in: right
         ! shape, already mapped, read-only in the callee, never indexed
         ! when `use_por` is .false.
         renorm_skip_walls = (bc_w_tag == OBC_WALL .and. bc_e_tag == OBC_WALL)
         ! FOUR branches, one per (porous, closed-faces) combination: each
         ! 3-D weight has to reach the callee as a full-size, device-present
         ! explicit-shape actual, and the knob-off ones hand over a read-only
         ! PPM edge buffer as the inert stand-in (see the `por` / `open_f`
         ! docstrings).  The DEFAULT path is the last branch and is textually
         ! the call this routine has always made, so it stays bit-identical.
         if (metrics%use_porous) then
            if (metrics%use_closed_faces) then
               call renormalise_zonal_flux_to_uhbt(grid, metrics, this, ms, uhbt, dt, &
                                                   skip_walls=renorm_skip_walls, &
                                                   has_west=has_w_flux, has_east=has_e_flux, &
                                                   visc_rem=visc_rem, u_cor=u_cor, &
                                                   use_por=.true., por=metrics%por_face_area_u, &
                                                   use_open=.true., open_f=metrics%open_u)
            else
               call renormalise_zonal_flux_to_uhbt(grid, metrics, this, ms, uhbt, dt, &
                                                   skip_walls=renorm_skip_walls, &
                                                   has_west=has_w_flux, has_east=has_e_flux, &
                                                   visc_rem=visc_rem, u_cor=u_cor, &
                                                   use_por=.true., por=metrics%por_face_area_u, &
                                                   use_open=.false., open_f=this%h_face_right_x%data)
            end if
         else if (metrics%use_closed_faces) then
            call renormalise_zonal_flux_to_uhbt(grid, metrics, this, ms, uhbt, dt, &
                                                skip_walls=renorm_skip_walls, &
                                                has_west=has_w_flux, has_east=has_e_flux, &
                                                visc_rem=visc_rem, u_cor=u_cor, &
                                                use_por=.false., por=this%h_face_left_x%data, &
                                                use_open=.true., open_f=metrics%open_u)
         else
            call renormalise_zonal_flux_to_uhbt(grid, metrics, this, ms, uhbt, dt, &
                                                skip_walls=renorm_skip_walls, &
                                                has_west=has_w_flux, has_east=has_e_flux, &
                                                visc_rem=visc_rem, u_cor=u_cor, &
                                                use_por=.false., por=this%h_face_left_x%data, &
                                                use_open=.false., open_f=this%h_face_right_x%data)
         end if
      end if
   end subroutine continuity_zonal_flux