continuity_tracer_step_split Subroutine

   subroutine continuity_tracer_step_split(grid, metrics, this, ms, dt, uhbt, vhbt, bc, mle, mle_fold_active, &
                                           tracer_mode, h_min, visc_rem_u, visc_rem_v, u_cor, v_cor)
      !! Production entry point for the directionally-split
      !! continuity + tracer step.  Interleaves the two so the
      !! CWC discrete theorem holds in the split form:
      !!
      !!   1. zonal_flux        — Φx from h^n
      !!   2. tracer_advect_zonal — hTr ← hTr - dt·∂(Φx·T)/∂x at h^n
      !!   3. apply_zonal       — h ← h^n - dt·∂Φx/∂x  (= h^*)
      !!   4. meridional_flux   — Φy from h^*
      !!   5. tracer_advect_meridional — hTr ← hTr - dt·∂(Φy·T)/∂y at h^*
      !!   6. apply_meridional  — h ← h^* - dt·∂Φy/∂y  (= h^{n+1})
      !!
      !! Uniform T preserved: after step 2, hTr = (h - dt·div_x)·T;
      !! after step 3, h = h - dt·div_x, so hTr/h = T still.  After
      !! step 5, hTr = (h^* - dt·div_y)·T = h^{n+1}·T.  After step 6,
      !! hTr/h = T.  Same CWC theorem as the unsplit form, lifted
      !! per direction.
      !!
      !! `flux_h_layer` ends the step holding the total horizontal
      !! divergence (sum of x and y substeps) — that's what the
      !! vertical-advection kernel consumes for w_interface.
      !!
      !! Optional `uhbt, vhbt`: time-mean barotropic-substep transports.  When
      !! supplied, the per-layer mass fluxes are renormalised so
      !! `Σ_k Φx_k = uhbt` and `Σ_k Φy_k = vhbt`, making the slow
      !! continuity advance `h_layer` consistently with the fast
      !! loop's `η_end` — MOM6's split-explicit pattern.  The same
      !! constrained fluxes feed tracer advection, so per-column
      !! `T = hTr/h` stays uniform under the constraint.
      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
      real(wp), intent(in), optional :: uhbt(:, :)
      real(wp), intent(in), optional :: vhbt(:, :)
      type(ocean_bc_state_t), intent(in), optional :: bc
         !! When present, per-edge OBC tags gate the wall-zero step
         !! inside the flux kernels.  OBC_WALL keeps the Phase 3
         !! closure; OBC_OPEN (and other non-wall tags) leaves the
         !! computed mass flux at the wall face for the downstream
         !! transport.  Absent ⇒ closed-wall everywhere.
      type(ocean_mle_t), intent(in), optional :: mle
         !! Fox-Kemper mixed-layer-eddy transports (B5).  When present
         !! and enabled, `mle%uhml`/`vhml` are folded into the per-layer
         !! mass fluxes AFTER each direction's flux fill and BEFORE the
         !! matching tracer advect + divergence — so the augmented flux
         !! transports both h and tracers (conservative; velocity
         !! untouched).  Absent / disabled ⇒ bit-identical no-op.
      logical, intent(in), optional :: mle_fold_active
         !! Gates the Fox-Kemper fold to the THERMO cadence.  Absent or
         !! `.true.` ⇒ the fold applies (bit-identical default — the case
         !! at `dt_therm_ratio = 1`, where every step is a thermo step).
         !! `.false.` skips the fold so the stale FK transports (computed
         !! once per thermo interval) are NOT re-applied on the
         !! intervening non-thermo outer steps when `dt_therm_ratio > 1`.
      ! assumed-shape-ok: pure passthroughs to the renormaliser.
      real(wp), intent(in), optional :: visc_rem_u(:, :, :), visc_rem_v(:, :, :)
         !! Per-layer viscous remnant gamma_k on east / north faces.  Forwarded
         !! to the flux renormalisers, where it weights the barotropic
         !! increment (MOM6 `u_cor = u + du*visc_rem`).  Absent => gamma == 1,
         !! bit-identical.
      ! assumed-shape-ok: pure passthroughs.
      real(wp), intent(inout), optional :: u_cor(:, :, :), v_cor(:, :, :)
         !! MOM6 `u_cor`/`v_cor` destinations — the step TIME-MEAN velocity
         !! (`u_av`/`v_av`), never the prognostic.  Absent => flux-only.
      integer, intent(in), optional :: tracer_mode
         !! Phase 2 (6b) windowed-advection mode.  `TR_MODE_ADVECT`
         !! (default, absent) ⇒ the historical fused path: advance h AND
         !! advect tracers each call (bit-identical to pre-6b).
         !! `TR_MODE_ACCUMULATE` ⇒ advance h, accumulate
         !! `0.5·mass_flux·dt` into `this%uhtr/vhtr` (one += per RK2
         !! stage, weight 0.5 baked in — closes the reconstruction against
         !! the RK2-averaged h), and SKIP the per-step tracer advect so
         !! `hTr` stays frozen until the boundary drain.
      real(wp), intent(in), optional :: h_min
         !! Phase-1 Lagrangian minimum-thickness floor (m). When > 0, passed
         !! to `continuity_apply_zonal`/`_meridional` to clamp h_new >= h_min.
         !! Absent or 0 ⇒ off ⇒ bit-identical.

      integer :: it
      logical :: per_x, per_y, do_mle_fold, fold_wall
      integer :: nx, ny, nz, nx_phys, ny_phys, nghost
      integer :: mode
      integer :: ii, jj, kk
      integer :: it_cw
      integer :: bc_w_tag, bc_e_tag, bc_s_tag, bc_n_tag
      real(wp) :: h_min_use

      mode = TR_MODE_ADVECT
      if (present(tracer_mode)) mode = tracer_mode

      h_min_use = 0.0_wp
      if (present(h_min)) h_min_use = h_min

      ! P2 positive-definite limiter: reset the per-call limited-face counter
      ! (accumulated across the zonal + meridional passes below).  Host scalar.
      this%n_limited_step = 0

      per_x = .false.
      per_y = .false.
      if (present(bc)) then
         per_x = bc%periodic_x .and. .not. ocean_halo_is_decomposed_x()
         per_y = bc%periodic_y .and. .not. ocean_halo_is_decomposed_y()
      end if
      ! Fox-Kemper fold defaults ON (bit-identical for callers that do not
      ! pass the gate); the dyn step passes `is_thermo_step()` to suppress
      ! the fold on non-thermo steps when dt_therm_ratio > 1.
      do_mle_fold = .true.
      if (present(mle_fold_active)) do_mle_fold = mle_fold_active
      ! Whether the MLE bolus fold actually contributes this call.  When it
      ! does, the augmented flux must be re-closed at no-normal-flow WALL
      ! faces (the fold adds bolus transport at every face, including the
      ! physical wall the resolved flux already zeroed — otherwise the bolus
      ! bleeds tracer mass into the ghost halo across the wall).  GM is no
      ! longer folded: it is its own operator (`continuity_gm_apply`).
      fold_wall = .false.
      if (present(mle)) fold_wall = fold_wall .or. mle%enable
      fold_wall = fold_wall .and. do_mle_fold
      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nx_phys = grid%nx_phys
      ny_phys = grid%ny_phys
      nghost = grid%nghost

      ! Windowed-advect concentration hold (step 1 of 2).  In
      ! TR_MODE_ACCUMULATE the horizontal tracer advect is deferred to the
      ! end-of-window drain, so `hTr` must not move — but `h_layer` does,
      ! every stage.  Snapshot the pre-continuity thickness so the paired
      ! rescale at the bottom of this routine can hold `T = hTr/h_layer`
      ! fixed instead of holding `hTr` fixed.  `hprev_work` is idle here:
      ! the drain is the only other consumer and it runs at outer-step end.
      if (mode == TR_MODE_ACCUMULATE) then
         ! First accumulate stage of a window: latch the window-start
         ! thickness the drain will use to undo the hold exactly.
         if (.not. this%hTr_holds_conc) then
            call drain_copy_3d(nx, ny, nz, ms%h_layer, this%h_win_start)
         end if
         call drain_copy_3d(nx, ny, nz, ms%h_layer, this%hprev_work)
      end if

      if (present(uhbt) .and. present(bc)) then
         call continuity_zonal_flux(grid, metrics, this, ms, dt, uhbt=uhbt, bc=bc, &
                                    visc_rem=visc_rem_u, u_cor=u_cor)
      else if (present(uhbt)) then
         call continuity_zonal_flux(grid, metrics, this, ms, dt, uhbt=uhbt, &
                                    visc_rem=visc_rem_u, u_cor=u_cor)
      else if (present(bc)) then
         call continuity_zonal_flux(grid, metrics, this, ms, dt, bc=bc)
      else
         call continuity_zonal_flux(grid, metrics, this, ms, dt)
      end if
      ! FK MLE fold (B5): add uhml into the zonal mass flux so it advects
      ! both h and tracers and enters the divergence.  No-op if absent.
      ! NOTE: this fold runs every dynamics call; `mle_compute_transports`
      ! runs only at thermo cadence.  When dt_therm_ratio > 1 the stale
      ! uhml/vhml are re-folded on non-thermo steps (over-applies FK at
      ! dynamics cadence — known limitation, safe/conservative via
      ! sum_k a(k) = 0; to be gated when sub-thermo cadence is exercised).
      if (present(mle) .and. do_mle_fold) then
         call mle_fold_x(mle, ms%mass_flux_x_layer, nx + 1, ny, nz)
      end if
      ! No-normal-flow wall closure for the bolus fold (mirrors the
      ! resolved-flux wall zeroing in `continuity_zonal_flux`): the MLE
      ! fold above adds `uhml` at the physical WALL faces, which the
      ! resolved flux had zeroed.  Without re-zeroing, the bolus transports
      ! tracer mass across the wall into the ghost halo and the
      ! physical-domain `sum(hTr)` drifts.  BC-aware: periodic / open edges
      ! keep the folded transport.
      if (fold_wall) then
         ! Re-close the physical walls after the MLE fold.  An MPI seam
         ! (`has_*` false) is not a wall: zeroing it cut every decomposed
         ! Fox-Kemper run's transport at the rank seams (the tile edge
         ! is an interior face the neighbour computes identically).
         bc_w_tag = OBC_WALL
         bc_e_tag = OBC_WALL
         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)
            if (.not. bc%has_west) bc_w_tag = OBC_PERIODIC
            if (.not. bc%has_east) bc_e_tag = OBC_PERIODIC
         end if
         do concurrent(kk=1:nz, jj=1:ny)
            if (bc_w_tag == OBC_WALL) ms%mass_flux_x_layer(nghost + 1, jj, kk) = 0.0_wp
            if (bc_e_tag == OBC_WALL) ms%mass_flux_x_layer(nghost + nx_phys + 1, jj, kk) = 0.0_wp
         end do
      end if
      ! P2 positive-definite outflux limiter (zonal): scale the OUTGOING
      ! east-face mass fluxes so no donor drains below h_lim.  Applied to the
      ! FOLDED total (after the MLE bolus fold + wall closure) so the
      ! h-apply, tracer advect, uhtr accumulation, and the corrector's
      ! `use_state_fluxes` reads all consume the SAME limited flux (D3).
      ! Off ⇒ skipped ⇒ bit-identical.
      if (this%positive_definite) then
         ! v1.1: forward u_cor so the limiter re-scales the captured
         ! transport-matched velocity (the u_av family) by the same θ as
         ! the flux — optional-forwarding propagates absence.  Keeps u_av
         ! consistent with the LIMITED fluxes the pred_corr corrector's
         ! use_state_fluxes CorAdCalc transports with (the flux↔velocity
         ! match is load-bearing).
         call pd_limit_zonal_impl(nx, ny, nz, dt, this%h_lim, metrics%iareaT, &
                                  ms%h_layer, ms%mass_flux_x_layer, &
                                  this%pd_theta%data, this%n_limited_step, &
                                  u_cor=u_cor)
      end if
      if (mode == TR_MODE_ACCUMULATE) then
         ! Windowed mode: accumulate this stage's zonal mass flux (with
         ! the RK2 0.5 weight baked in) and SKIP the tracer advect so
         ! hTr stays frozen.  hprev = h^{n+1} + div(uhtr) then closes
         ! the reconstruction against the RK2-averaged h.
         call accumulate_flux_x(nx + 1, ny, nz, dt, ms%mass_flux_x_layer, this%uhtr)
      else if (mode /= TR_MODE_NONE .and. present(bc)) then
         call tracer_advect_zonal(grid, metrics, this, ms, dt, bc=bc)
      else if (mode /= TR_MODE_NONE) then
         call tracer_advect_zonal(grid, metrics, this, ms, dt)
      end if
      if (h_min_use > 0.0_wp) then
         call continuity_apply_zonal(grid, metrics, ms, dt, h_min=h_min_use)
      else
         call continuity_apply_zonal(grid, metrics, ms, dt)
      end if

      ! Mid-Lie-split MPI seam exchange (O3 correctness fix): after the zonal
      ! apply updates h_layer (and hTr in ADVECT mode) the meridional flux
      ! reconstruction reads h ghost columns that the NEIGHBOUR rank's zonal
      ! apply has updated — but those ghosts were never re-exchanged.  Without
      ! this exchange an x-rank seam leaks ~3e-9/day global mass.
      ! D0-unconditional: single-rank non-periodic => no-op, single-rank
      ! periodic => local wrap, multi-rank => messages.  Must run BEFORE the
      ! periodic wrap below so both exchanges see the same post-zonal state.
      !
      ! NOTE: this exchange is also inside the "ocean_continuity" compute region
      ! opened by the dyn caller (rdb_ocean_dyn.F90); ocean_comms_ml here isolates
      ! the comm share, so sums of compute+comms slightly over-close by this term.
      ! This is a known, accepted double-attribution — no stop/restart of the
      ! outer region from inside this module.
      call profiler_start("ocean_comms_ml")
      call ocean_halo_centre(ms%h_layer, nz)
      if (allocated(ms%tracers)) then
         do it = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it)%hTr)) cycle
            call ocean_halo_centre(ms%tracers(it)%hTr, nz)
         end do
      end if
      call profiler_stop("ocean_comms_ml")

      ! Mid-Lie-split ghost wrap (design §1.5): re-wrap h_layer (and
      ! tracers) after the zonal apply so the meridional reconstruction
      ! reads current ghost values.  Both per_x and per_y checked: even
      ! for periodic-x only, per_y ghosts can inherit stale values
      ! accumulated under the zonal update.  Two cheap DC kernels, no
      ! correctness traps.
      if (per_x .or. per_y) then
         ! Batched async wrap (queue 1): h_layer + every tracer issued without
         ! per-call sync, then synced ONCE below — pipelines the tiny ghost-slab
         ! launches (otherwise launch-latency-bound).  Wait before the fold,
         ! which reads these wrapped ghosts.
         call ocean_periodic_wrap_centre_3d(ms%h_layer, nx, ny, nz, &
                                            nx_phys, ny_phys, nghost, per_x, per_y, no_wait=.true.)
         if (allocated(ms%tracers)) then
            do it = 1, size(ms%tracers)
               if (.not. allocated(ms%tracers(it)%hTr)) cycle
               call ocean_periodic_wrap_centre_3d(ms%tracers(it)%hTr, nx, ny, nz, &
                                                  nx_phys, ny_phys, nghost, per_x, per_y, no_wait=.true.)
            end do
         end if
         !$acc wait(1)
      end if
      ! Mid-Lie-split north fold (Appendix A): re-fold the centre fields
      ! (h_layer + tracers) AFTER the periodic wrap so the meridional
      ! reconstruction reads fold-consistent north ghosts.  No-op when not
      ! folding.
      if (present(bc)) call ocean_fold_wrap_centre_3d_state(grid, bc, ms)

      if (present(vhbt) .and. present(bc)) then
         call continuity_meridional_flux(grid, metrics, this, ms, dt, vhbt=vhbt, bc=bc, &
                                         visc_rem=visc_rem_v, v_cor=v_cor)
      else if (present(vhbt)) then
         call continuity_meridional_flux(grid, metrics, this, ms, dt, vhbt=vhbt, &
                                         visc_rem=visc_rem_v, v_cor=v_cor)
      else if (present(bc)) then
         call continuity_meridional_flux(grid, metrics, this, ms, dt, bc=bc)
      else
         call continuity_meridional_flux(grid, metrics, this, ms, dt)
      end if
      ! FK MLE fold (B5): add vhml into the meridional mass flux.
      ! Same thermo-cadence / dynamics-fold limitation as the zonal fold above.
      if (present(mle) .and. do_mle_fold) then
         call mle_fold_y(mle, ms%mass_flux_y_layer, nx, ny + 1, nz)
      end if
      ! No-normal-flow wall closure for the meridional bolus fold (see the
      ! zonal block above for the rationale).
      if (fold_wall) then
         ! MPI seams are not walls (see the zonal twin above).
         bc_s_tag = OBC_WALL
         bc_n_tag = OBC_WALL
         if (present(bc)) then
            bc_s_tag = ocean_bc_outer_face_tag(bc%south%bc_type)
            bc_n_tag = ocean_bc_outer_face_tag(bc%north%bc_type)
            if (.not. bc%has_south) bc_s_tag = OBC_PERIODIC
            if (.not. bc%has_north) bc_n_tag = OBC_PERIODIC
         end if
         do concurrent(kk=1:nz, ii=1:nx)
            if (bc_s_tag == OBC_WALL) ms%mass_flux_y_layer(ii, nghost + 1, kk) = 0.0_wp
            if (bc_n_tag == OBC_WALL) ms%mass_flux_y_layer(ii, nghost + ny_phys + 1, kk) = 0.0_wp
         end do
      end if
      ! P2 positive-definite outflux limiter (meridional): mirror of the zonal
      ! pass, on the post-zonal-apply h* availability.  Same D3 single-source
      ! scaling of the folded total.  Off ⇒ skipped ⇒ bit-identical.
      if (this%positive_definite) then
         ! v1.1: forward v_cor — see the zonal twin.
         call pd_limit_meridional_impl(nx, ny, nz, dt, this%h_lim, metrics%iareaT, &
                                       ms%h_layer, ms%mass_flux_y_layer, &
                                       this%pd_theta%data, this%n_limited_step, &
                                       v_cor=v_cor)
      end if
      ! Tripolar fold-line flux projection.  The fold-line row (north face
      ! of the last T-row, `rdb_ocean_fold` header) stores ONE physical face
      ! twice; its two flux slots were computed independently (PPM
      ! reconstruction, vhbt renormalisation, MLE bolus, limiter) and
      ! agree only up to rounding.  Project the FINAL flux antisymmetric
      ! (and refill the rows above it) before it touches h / hTr / the
      ! accumulated transport, so the cross-fold exchange telescopes: what
      ! leaves cell (i,nj) through its north face is exactly what enters
      ! cell (ni+1-i,nj).  No-op when not folding.  px > 1: no halo
      ! precedes this point, so the owner-routed exchange is what makes it
      ! exact (every value comes from the rank that owns the mirror face).
      if (present(bc)) then
         if (bc%north_fold) call ocean_fold_north_v_face(ms%mass_flux_y_layer, nx, ny + 1, nz, &
                                                         nx_phys, ny_phys, nghost)
      end if
      if (mode == TR_MODE_ACCUMULATE) then
         call accumulate_flux_y(nx, ny + 1, nz, dt, ms%mass_flux_y_layer, this%vhtr)
      else if (mode /= TR_MODE_NONE .and. present(bc)) then
         call tracer_advect_meridional(grid, metrics, this, ms, dt, bc=bc)
      else if (mode /= TR_MODE_NONE) then
         call tracer_advect_meridional(grid, metrics, this, ms, dt)
      end if
      if (h_min_use > 0.0_wp) then
         call continuity_apply_meridional(grid, metrics, ms, dt, h_min=h_min_use)
      else
         call continuity_apply_meridional(grid, metrics, ms, dt)
      end if
      ! Windowed-advect concentration hold (step 2 of 2).  Both applies have
      ! advanced `h_layer`; re-weight the frozen tracer content onto it so the
      ! concentration every downstream consumer reads (`ocean_eos_compute` →
      ! ρ → PGF, vdiff, hdiff, vertical advect, diagnostics) is EXACTLY the
      ! window-start value, as it is in MOM6 (whose prognostic `Tr%t` is a
      ! concentration and is therefore thickness-invariant for free).
      ! Without this, `T` drifts by the full window thickness divergence and
      ! the resulting grid-scale buoyancy error closes an exponentially
      ! growing EOS→PGF→divergence loop.  Composes with `rk2_average`:
      ! `hTr0 = T·h^n` and `hTr = T·h^(2)` average to `T·h^(n+1)`, so `T` is
      ! still exactly `T`.  `continuity_tracer_drain` converts back before it
      ! spends the accumulated transports.
      if (mode == TR_MODE_ACCUMULATE .and. allocated(ms%tracers)) then
         do it_cw = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it_cw)%hTr)) cycle
            if (.not. ms%tracers(it_cw)%do_horizontal_advection) cycle
            ! Budget dispatch (mirrors `tracer_advect_zonal`): heat/salt get
            ! the hold's own content change recorded so the console budget
            ! closes at EVERY report, not only on window boundaries.  Tracers
            ! without a budget slot take the budget-free twin.
            select case (ms%tracers(it_cw)%budget_id)
            case (TRACER_BUDGET_HEAT)
               call drain_rescale_hTr_budget(nx, ny, nz, ms%h_layer, this%hprev_work, &
                                             DRAIN_BUDGET_IN_STAGE_WEIGHT, &
                                             ms%tracers(it_cw)%hTr, ms%heat_budget_horiz_adv)
            case (TRACER_BUDGET_SALT)
               call drain_rescale_hTr_budget(nx, ny, nz, ms%h_layer, this%hprev_work, &
                                             DRAIN_BUDGET_IN_STAGE_WEIGHT, &
                                             ms%tracers(it_cw)%hTr, ms%salt_budget_horiz_adv)
            case default
               call drain_rescale_hTr(nx, ny, nz, ms%h_layer, this%hprev_work, &
                                      ms%tracers(it_cw)%hTr)
            end select
         end do
         this%hTr_holds_conc = .true.
      end if

      ! P3: fold this call's limited-face count into the running total (one
      ! add per split call, mirroring dyn%ntrunc_total).  n_limited_step is 0
      ! when positive_definite is off, so this is a no-op there.
      this%n_limited_total = this%n_limited_total + this%n_limited_step
   end subroutine continuity_tracer_step_split