redi_apply_flux_impl Subroutine

   subroutine redi_apply_flux_impl(nx, ny, nz, ns, dt, &
                                   nghost, nxp, nyp, wall_w, wall_e, wall_s, wall_n, &
                                   khtr_u, khtr_v, dy_cu, dx_cv, &
                                   idxCu, idyCv, areaT, h_layer, hTr_in, hTr, &
                                   uPoL, uPoR, uKoL, uKoR, uhEff, &
                                   vPoL, vPoR, vKoL, vKoR, vhEff, &
                                   uKb, uKt, vKb, vKt, use_open)
      !! Flat-impl Phase-B kernel for ONE tracer.  Cell-centric double-visit
      !! (no-scatter rule on the C-grid): cell (i,j) recomputes the
      !! along-neutral flux on each of its four bounding faces and accumulates
      !! ONLY into its own `dTr`; interior-face fluxes are computed twice but
      !! no thread writes a neighbour ⇒ race-free.
      !!   Face sign: the LEFT (west/south) cell of a face gets +Flx into
      !!   native layer nz+1-KoL; the RIGHT (east/north) cell gets -Flx into
      !!   nz+1-KoR.  The top-down→native flip k=nz+1-Ko is the only k-flip.
      !!   Flux = dT_layer * hEff * Coef, Coef_u = dt*khtr_u*dy_cu*idxCu;
      !!   divergence hTr(k) += dTr(k)/areaT (conservative).
      integer, intent(in) :: nx, ny, nz, ns
      real(wp), intent(in) :: dt
      integer, intent(in) :: nghost, nxp, nyp
      logical, intent(in) :: wall_w, wall_e, wall_s, wall_n
      real(wp), intent(in) :: khtr_u(nx + 1, ny)
      real(wp), intent(in) :: khtr_v(nx, ny + 1)
      real(wp), intent(in) :: dy_cu(nx + 1, ny)
      real(wp), intent(in) :: dx_cv(nx, ny + 1)
      real(wp), intent(in) :: idxCu(nx + 1, ny)
      real(wp), intent(in) :: idyCv(nx, ny + 1)
      real(wp), intent(in) :: areaT(nx, ny)
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: hTr_in(nx, ny, nz)
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(in) :: uPoL(nx + 1, ny, ns), uPoR(nx + 1, ny, ns)
      integer, intent(in) :: uKoL(nx + 1, ny, ns), uKoR(nx + 1, ny, ns)
      real(wp), intent(in) :: uhEff(nx + 1, ny, ns - 1)
      real(wp), intent(in) :: vPoL(nx, ny + 1, ns), vPoR(nx, ny + 1, ns)
      integer, intent(in) :: vKoL(nx, ny + 1, ns), vKoR(nx, ny + 1, ns)
      real(wp), intent(in) :: vhEff(nx, ny + 1, ns - 1)
      integer, intent(in) :: uKb(nx + 1, ny), uKt(nx + 1, ny)
         !! Phase-A u-face open windows (`1..nz` off the closed-face path).
      integer, intent(in) :: vKb(nx, ny + 1), vKt(nx, ny + 1)
         !! Phase-A v-face open windows.
      logical, intent(in) :: use_open
         !! z-level closed faces active: re-check window liveness.

      integer :: i, j, k
      integer :: wuf_w, wuf_e, wvf_s, wvf_n
      real(wp) :: dTr(NZ_STACK_MAX)
      real(wp) :: iaij

      ! No-normal-flow WALL faces of the physical domain.  The Redi flux must
      ! NOT cross these — else it bleeds tracer into the ghost halo.  A wall
      ! face is skipped by both adjacent cells, so the +flx/-flx pair never forms.
      wuf_w = nghost + 1
      wuf_e = nghost + nxp + 1
      wvf_s = nghost + 1
      wvf_n = nghost + nyp + 1

      ! Per-cell double-visit: each interior face is computed by both adjacent
      ! cells (race-free — every cell writes only its own dTr).  The heavy
      ! column locals live inside `redi_face_flux`, keeping this loop body's
      ! footprint to `dTr`/`iaij` (avoids the gfortran-15.1 dc-local corruption).
      do concurrent(j=1:ny, i=1:nx) local(k, dTr, iaij)
         do k = 1, nz
            dTr(k) = 0.0_wp
         end do
         ! WEST u-face (face i): this cell is the RIGHT column.
         if (i >= 2 .and. .not. ((wall_w .and. i == wuf_w) .or. (wall_e .and. i == wuf_e))) then
            call redi_face_flux(nz, ns, nx, ny, nx + 1, ny, h_layer, hTr_in, &
                                i - 1, j, i, j, i, j, uPoL, uPoR, uKoL, uKoR, uhEff, &
                                uKb(i, j), uKt(i, j), use_open, &
                                dt*khtr_u(i, j)*dy_cu(i, j)*idxCu(i, j), .false., dTr)
         end if
         ! EAST u-face (face i+1): this cell is the LEFT column.
         if (i <= nx - 1 .and. .not. ((wall_w .and. i + 1 == wuf_w) .or. (wall_e .and. i + 1 == wuf_e))) then
            call redi_face_flux(nz, ns, nx, ny, nx + 1, ny, h_layer, hTr_in, &
                                i, j, i + 1, j, i + 1, j, uPoL, uPoR, uKoL, uKoR, uhEff, &
                                uKb(i + 1, j), uKt(i + 1, j), use_open, &
                                dt*khtr_u(i + 1, j)*dy_cu(i + 1, j)*idxCu(i + 1, j), .true., dTr)
         end if
         ! SOUTH v-face (face j): this cell is the RIGHT (north) column.
         if (j >= 2 .and. .not. ((wall_s .and. j == wvf_s) .or. (wall_n .and. j == wvf_n))) then
            call redi_face_flux(nz, ns, nx, ny, nx, ny + 1, h_layer, hTr_in, &
                                i, j - 1, i, j, i, j, vPoL, vPoR, vKoL, vKoR, vhEff, &
                                vKb(i, j), vKt(i, j), use_open, &
                                dt*khtr_v(i, j)*dx_cv(i, j)*idyCv(i, j), .false., dTr)
         end if
         ! NORTH v-face (face j+1): this cell is the LEFT (south) column.
         if (j <= ny - 1 .and. .not. ((wall_s .and. j + 1 == wvf_s) .or. (wall_n .and. j + 1 == wvf_n))) then
            call redi_face_flux(nz, ns, nx, ny, nx, ny + 1, h_layer, hTr_in, &
                                i, j, i, j + 1, i, j + 1, vPoL, vPoR, vKoL, vKoR, vhEff, &
                                vKb(i, j + 1), vKt(i, j + 1), use_open, &
                                dt*khtr_v(i, j + 1)*dx_cv(i, j + 1)*idyCv(i, j + 1), .true., dTr)
         end if
         iaij = 1.0_wp/areaT(i, j)
         do k = 1, nz
            hTr(i, j, k) = hTr(i, j, k) + dTr(k)*iaij
         end do
      end do
   end subroutine redi_apply_flux_impl