redi_face_flux Subroutine

   pure subroutine redi_face_flux(nz, ns, nxc, nyc, nfa, nfb, h_layer, hTr_in, &
                                  iL, jL, iR, jR, fa, fb, &
                                  PoL, PoR, KoL, KoR, hEff, kb, kt, use_open, &
                                  coef, is_left, dTr)
      !$acc routine seq
      !! Accumulate ONE C-grid face's neutral-surface tracer flux into the
      !! owning cell's `dTr`.  Builds the left/right tracer columns from the
      !! READ-ONLY snapshot `hTr_in` (the live `hTr` is also written by the
      !! loop — reading it would be a do-concurrent read-write race) and loops
      !! the `ns-1` neutral sublayers.  The per-face column locals live in this
      !! frame to keep the caller's loop-body footprint small.
      !! `is_left`: this cell is the LEFT (west/south) column ⇒ `+flx` into
      !! native layer `nz+1-KoL`; else the RIGHT column ⇒ `-flx` into
      !! `nz+1-KoR`.  `(iL,jL)`/`(iR,jR)` index the columns; `(fa,fb)` the faces.
      !! `kb..kt` is the face's Phase-A open window (`1..nz` off the z-level
      !! closed-face path): the tracer columns are reconstructed on it alone
      !! and the full-frame `Ko` are read in the window frame (`Ko - nz +
      !! kt`).  An empty window, or (`use_open`) one whose layers are no
      !! longer all live on both sides, contributes nothing — both cells of
      !! the face take the same decision from the same `h`.
      integer, intent(in) :: nz, ns, nxc, nyc, nfa, nfb
      integer, intent(in) :: iL, jL, iR, jR, fa, fb
      real(wp), intent(in) :: h_layer(nxc, nyc, nz), hTr_in(nxc, nyc, nz)
      real(wp), intent(in) :: PoL(nfa, nfb, ns), PoR(nfa, nfb, ns)
      integer, intent(in) :: KoL(nfa, nfb, ns), KoR(nfa, nfb, ns)
      real(wp), intent(in) :: hEff(nfa, nfb, ns - 1)
      integer, intent(in) :: kb, kt
      logical, intent(in) :: use_open
      real(wp), intent(in) :: coef
      logical, intent(in) :: is_left
      real(wp), intent(inout) :: dTr(nz)

      integer :: k, ks, knat, nk, koff
      real(wp) :: hcL(NZ_STACK_MAX), trcL(NZ_STACK_MAX)
      real(wp) :: hcR(NZ_STACK_MAX), trcR(NZ_STACK_MAX)
      real(wp) :: TlL(NZ_STACK_MAX), TiL(NZ_STACK_MAX + 1), aLL(NZ_STACK_MAX), aRL(NZ_STACK_MAX)
      real(wp) :: TlR(NZ_STACK_MAX), TiR(NZ_STACK_MAX + 1), aLR(NZ_STACK_MAX), aRR(NZ_STACK_MAX)
      real(wp) :: dtdiff, flx

      if (kt < kb) return
      if (use_open) then
         do k = kb, kt
            if (.not. (rdb_vl_is_live(h_layer(iL, jL, k)) .and. &
                       rdb_vl_is_live(h_layer(iR, jR, k)))) return
         end do
      end if
      nk = kt - kb + 1
      koff = nz - kt
      do k = 1, nz
         hcL(k) = h_layer(iL, jL, k)
         trcL(k) = hTr_in(iL, jL, k)
         hcR(k) = h_layer(iR, jR, k)
         trcR(k) = hTr_in(iR, jR, k)
      end do
      call redi_tracer_column(kb, kt, hcL, trcL, TlL, TiL, aLL, aRL)
      call redi_tracer_column(kb, kt, hcR, trcR, TlR, TiR, aLR, aRR)
      do ks = 1, ns - 1
         if (hEff(fa, fb, ks) /= 0.0_wp) then
            dtdiff = redi_sublayer_dT(nk, KoL(fa, fb, ks) - koff, KoL(fa, fb, ks + 1) - koff, &
                                      KoR(fa, fb, ks) - koff, KoR(fa, fb, ks + 1) - koff, &
                                      PoL(fa, fb, ks), PoL(fa, fb, ks + 1), &
                                      PoR(fa, fb, ks), PoR(fa, fb, ks + 1), &
                                      TlL, TiL, aLL, aRL, TlR, TiR, aLR, aRR)
            flx = dtdiff*hEff(fa, fb, ks)*coef
            if (is_left) then
               knat = nz + 1 - KoL(fa, fb, ks)
               dTr(knat) = dTr(knat) + flx
            else
               knat = nz + 1 - KoR(fa, fb, ks)
               dTr(knat) = dTr(knat) - flx
            end if
         end if
      end do
   end subroutine redi_face_flux