redi_apply_flux Subroutine

   subroutine redi_apply_flux(grid, metrics, this, ms, dt, khtr_u_ext, khtr_v_ext, bc)
      !! Public Phase-B entry: apply the neutral-diffusion tracer update for
      !! every registered tracer.  No-op if absent / uninit / disabled / zero
      !! diffusivity.  Run at thermo cadence after `redi_calc_coeffs` and the
      !! along-coordinate `tracer_hdiff` (Redi augments it).  Per-face KhTr
      !! comes from `khtr_u_ext`/`khtr_v_ext` (VarMix) when supplied, else the
      !! scalar `this%khtr` broadcast onto every face.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_redi_t), intent(inout) :: this
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      real(wp), intent(in), optional :: khtr_u_ext(:, :)
         !! VarMix per-face KhTr at u-faces `(nx+1, ny)` (m^2/s).
      real(wp), intent(in), optional :: khtr_v_ext(:, :)
         !! VarMix per-face KhTr at v-faces `(nx, ny+1)` (m^2/s).
      type(ocean_bc_state_t), intent(in), optional :: bc
         !! Per-edge OBC tags.  No along-isopycnal flux crosses a
         !! no-normal-flow (`OBC_WALL`) physical-domain boundary face — else
         !! Redi bleeds tracer into the ghost halo.  An OPEN (tracer-open)
         !! edge KEEPS its face flux, read against the OBC-filled ghost column
         !! — MOM6 `neutral_diffusion` gates its faces on `G%mask2dCu`, which
         !! `open_boundary_impose_land_mask` leaves at 1 on an open segment's
         !! normal face (it zeroes only `OBCmaskCu` there), so neutral
         !! diffusion exchanges tracer with the exterior; that exchange is
         !! booked in `*_budget_hdiff` below so the closed budget stays
         !! closed.  An MPI seam (`has_* = .false.`) is never a wall.
         !! Absent ⇒ all edges WALL.
      integer :: nx, ny, nz, it
      integer :: nghost, nxp, nyp
      logical :: use_ext, wall_w, wall_e, wall_s, wall_n

      if (.not. this%is_init) return
      if (.not. this%enable) return
      if (.not. allocated(ms%h_layer)) return
      if (.not. allocated(ms%tracers)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      if (this%nz_ml /= nz) return
      nghost = grid%nghost
      nxp = grid%nx_phys
      nyp = grid%ny_phys
      ! Which physical-domain edges are no-normal-flow walls (default closed).
      ! The edge tags are GLOBAL: under a domain decomposition a tile edge
      ! that is not a physical domain edge (`.not. bc%has_*`) is an MPI
      ! seam, which the neutral flux must cross exactly as an interior face
      ! does — keyed on the tag alone, a wall-tagged edge would close every
      ! seam parallel to it (`test_ocean_decomp_bitid_mpi`, 1x2 split).
      ! Single rank: `has_* = .true.` ⇒ unchanged.
      wall_w = .true.
      wall_e = .true.
      wall_s = .true.
      wall_n = .true.
      if (present(bc)) then
         wall_w = (ocean_bc_outer_face_tag(bc%west%bc_type) == OBC_WALL) .and. bc%has_west
         wall_e = (ocean_bc_outer_face_tag(bc%east%bc_type) == OBC_WALL) .and. bc%has_east
         wall_s = (ocean_bc_outer_face_tag(bc%south%bc_type) == OBC_WALL) .and. bc%has_south
         wall_n = (ocean_bc_outer_face_tag(bc%north%bc_type) == OBC_WALL) .and. bc%has_north
      end if

      use_ext = present(khtr_u_ext) .and. present(khtr_v_ext)
      ! No diffusivity anywhere => nothing to do (scalar path only; the
      ! VarMix field may be non-zero even when the scalar floor is 0).
      if (.not. use_ext .and. this%khtr <= 0.0_wp) return

      ! Resolve the per-face KhTr into the slot fields (device-resident).
      if (use_ext) then
         call redi_face_copy(nx + 1, ny, khtr_u_ext, this%khtr_u)
         call redi_face_copy(nx, ny + 1, khtr_v_ext, this%khtr_v)
      else
         call redi_face_const(nx + 1, ny, this%khtr, this%khtr_u)
         call redi_face_const(nx, ny + 1, this%khtr, this%khtr_v)
      end if

      do it = 1, size(ms%tracers)
         ! Snapshot this tracer's hTr so the double-visit gather reads the
         ! pre-step field (order-independent — see `tr_snap`).
         call redi_snapshot(nx, ny, nz, ms%tracers(it)%hTr, this%tr_snap)
         call redi_apply_flux_impl(nx, ny, nz, this%nsurf, dt, &
                                   nghost, nxp, nyp, wall_w, wall_e, wall_s, wall_n, &
                                   this%khtr_u, this%khtr_v, &
                                   metrics%dy_cu, metrics%dx_cv, &
                                   metrics%idxCu, metrics%idyCv, metrics%areaT, &
                                   ms%h_layer, this%tr_snap, ms%tracers(it)%hTr, &
                                   this%uPoL, this%uPoR, this%uKoL, this%uKoR, this%uhEff, &
                                   this%vPoL, this%vPoR, this%vKoL, this%vKoR, this%vhEff, &
                                   this%uKb, this%uKt, this%vKb, this%vKt, &
                                   metrics%use_closed_faces)
         ! Closed-budget bookkeeping (salt / heat): fold the realised Redi
         ! change into the same lateral-diffusion accumulator `tracer_hdiff`
         ! uses.  Interior faces cancel in the interior sum, so what survives
         ! is exactly the flux through the OPEN physical faces — the term the
         ! console `out` column must carry.  `hTr - tr_snap` is the realised
         ! increment and is exact (Sterbenz: the two are within a factor 2).
         select case (ms%tracers(it)%budget_id)
         case (TRACER_BUDGET_SALT)
            call redi_budget_accumulate(nx, ny, nz, ms%tracers(it)%hTr, this%tr_snap, &
                                        ms%salt_budget_hdiff)
         case (TRACER_BUDGET_HEAT)
            call redi_budget_accumulate(nx, ny, nz, ms%tracers(it)%hTr, this%tr_snap, &
                                        ms%heat_budget_hdiff)
         case default
            ! Passive tracers carry no closed budget.
         end select
      end do
   end subroutine redi_apply_flux