tracer_hdiff_one_impl Subroutine

   pure subroutine tracer_hdiff_one_impl(nx, ny, nz, nghost, nx_phys, ny_phys, dt, kappa, &
                                         dy_cu, dx_cv, idxCu, idyCv, iareaT, &
                                         h, hTr, T_centre, F_x_face, F_y_face, &
                                         wall_w, wall_e, wall_s, wall_n, &
                                         open_u, open_v, budget)
      !! Four-pass flat-impl horizontal-Laplacian tracer diffusion in
      !! conservative curvilinear form (design §2, mirrors continuity):
      !!
      !!   1. Compute T = hTr / h at every cell centre (vanishing
      !!      layer → T = 0).
      !!   2. East-face TRANSPORT flux [tracer·m³/s]:
      !!        F_x(i) = kappa · 0.5·(h(i-1)+h(i)) · (T(i)-T(i-1))
      !!                 · idxCu(i) · dy_cu(i)
      !!      PHYSICAL wall faces (i = nghost+1, i = nghost+nx_phys+1):
      !!      zeroed when `wall_w`/`wall_e` (mirrors
      !!      `continuity_zonal_flux`'s host-resolved has_*/OBC_WALL
      !!      gate).  ARRAY-bound faces (i = 1, i = nx+1) are always
      !!      zeroed too — the pass-4 divergence reads F_x_face(1,...)
      !!      and F_x_face(nx+1,...) at the domain edge cells even
      !!      when the physical wall sits inside the ghost band, so
      !!      those two faces must stay initialised.
      !!   3. North-face flux, mirror of pass 2 (idyCv · dx_cv).
      !!   4. Divergence into hTr:
      !!        hTr(i,j,k) += dt · [(F_x(i+1)-F_x(i))
      !!                            + (F_y(j+1)-F_y(j))] · iareaT(i,j)
      !!
      !! On uniform Cartesian `dy_cu=dy`, `idxCu=1/dx`,
      !! `iareaT=1/(dx·dy)`, so the form collapses to the old
      !! `Δ(kappa·h·ΔT/dx)/dx` to round-off.  Each pass writes a
      !! different buffer than it reads, so `do concurrent` is
      !! race-free.  Constancy: uniform T → zero face fluxes → hTr
      !! unchanged.  Conservation over the PHYSICAL domain: closed
      !! physical-wall + flux-form divergence over the area-weighted
      !! cells integrates to zero exactly (see `tracer_hdiff`'s
      !! host-side wall-flag resolution for the OBC/MPI-seam gate).
      !!
      !! Stability (per cell, explicit forward-Euler):
      !!   kappa · dt · (idxT² + idyT²) <= 0.5.
      integer, intent(in) :: nx, ny, nz, nghost, nx_phys, ny_phys
      real(wp), intent(in) :: dt, kappa
      real(wp), intent(in) :: dy_cu(nx + 1, ny), dx_cv(nx, ny + 1)
      real(wp), intent(in) :: idxCu(nx + 1, ny), idyCv(nx, ny + 1)
      real(wp), intent(in) :: iareaT(nx, ny)
      real(wp), intent(in) :: h(nx, ny, nz)
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: T_centre(nx, ny, nz)
      real(wp), intent(inout) :: F_x_face(nx + 1, ny, nz)
      real(wp), intent(inout) :: F_y_face(nx, ny + 1, nz)
      logical, intent(in) :: wall_w, wall_e, wall_s, wall_n
      real(wp), intent(in), optional :: open_u(nx + 1, ny, nz)
         !! Per-layer 0/1 u-face open mask
         !! (`&vcoord_nml zfixed_closed_faces`).  A CLOSED face is a
         !! z-level WALL for that layer, so it carries no diffusive
         !! tracer flux either — the same statement continuity makes
         !! about mass.
         !!
         !! ABSENT (the default path) ⇒ no masking pass is generated and
         !! every expression below is byte-identical to the un-masked
         !! form.  It is OPTIONAL rather than a `use_open` + inert
         !! stand-in pair precisely because the `(1,1,1)` placeholder
         !! must never reach an explicit-shape dummy, and this routine
         !! has no full-size read-only array of its own to lend (the
         !! `F_*_face` buffers it would otherwise borrow are its own
         !! `intent(inout)` scratch, so lending them would alias).
      real(wp), intent(in), optional :: open_v(nx, ny + 1, nz)
         !! v-face twin.  Present iff `open_u` is.
      real(wp), intent(inout), optional :: budget(nx, ny, nz)

      integer :: i, j, k
      real(wp) :: h_face, delta

      ! ---- Pass 1: T = hTr/h at centres ----
      do concurrent(k=1:nz, j=1:ny, i=1:nx)
         if (h(i, j, k) > 0.0_wp) then
            T_centre(i, j, k) = hTr(i, j, k)/h(i, j, k)
         else
            T_centre(i, j, k) = 0.0_wp
         end if
      end do

      ! ---- Pass 2: east-face transport flux ----
      do concurrent(k=1:nz, j=1:ny, i=2:nx) local(h_face)
         h_face = 0.5_wp*(h(i - 1, j, k) + h(i, j, k))
         F_x_face(i, j, k) = kappa*h_face* &
                             (T_centre(i, j, k) - T_centre(i - 1, j, k))* &
                             idxCu(i, j)*dy_cu(i, j)
      end do
      ! z-level closed faces: a separate host-gated pass so the loop above
      ! is textually unchanged with the knob off.
      if (present(open_u)) then
         do concurrent(k=1:nz, j=1:ny, i=2:nx)
            F_x_face(i, j, k) = F_x_face(i, j, k)*open_u(i, j, k)
         end do
      end if
      ! Array-bound faces: always zeroed (pass-4 divergence at the
      ! domain-edge cells reads them even when the physical wall is
      ! elsewhere inside the ghost band).
      do concurrent(k=1:nz, j=1:ny)
         F_x_face(1, j, k) = 0.0_wp
         F_x_face(nx + 1, j, k) = 0.0_wp
      end do
      ! Physical wall faces: zero only when the edge is actually a
      ! closed wall here (single-rank all-wall default; OBC/MPI-seam
      ! edges keep the computed flux read from a correctly-filled
      ! ghost column — see `tracer_hdiff`).
      do concurrent(k=1:nz, j=1:ny)
         if (wall_w) F_x_face(nghost + 1, j, k) = 0.0_wp
         if (wall_e) F_x_face(nghost + nx_phys + 1, j, k) = 0.0_wp
      end do

      ! ---- Pass 3: north-face transport flux ----
      do concurrent(k=1:nz, j=2:ny, i=1:nx) local(h_face)
         h_face = 0.5_wp*(h(i, j - 1, k) + h(i, j, k))
         F_y_face(i, j, k) = kappa*h_face* &
                             (T_centre(i, j, k) - T_centre(i, j - 1, k))* &
                             idyCv(i, j)*dx_cv(i, j)
      end do
      ! z-level closed faces: see the zonal twin.
      if (present(open_v)) then
         do concurrent(k=1:nz, j=2:ny, i=1:nx)
            F_y_face(i, j, k) = F_y_face(i, j, k)*open_v(i, j, k)
         end do
      end if
      do concurrent(k=1:nz, i=1:nx)
         F_y_face(i, 1, k) = 0.0_wp
         F_y_face(i, ny + 1, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, i=1:nx)
         if (wall_s) F_y_face(i, nghost + 1, k) = 0.0_wp
         if (wall_n) F_y_face(i, nghost + ny_phys + 1, k) = 0.0_wp
      end do

      ! ---- Pass 4: apply area-weighted divergence ----
      if (present(budget)) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx) local(delta)
            delta = dt*( &
                    (F_x_face(i + 1, j, k) - F_x_face(i, j, k)) + &
                    (F_y_face(i, j + 1, k) - F_y_face(i, j, k)))*iareaT(i, j)
            hTr(i, j, k) = hTr(i, j, k) + delta
            budget(i, j, k) = budget(i, j, k) + delta
         end do
      else
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            hTr(i, j, k) = hTr(i, j, k) + dt*( &
                           (F_x_face(i + 1, j, k) - F_x_face(i, j, k)) + &
                           (F_y_face(i, j + 1, k) - F_y_face(i, j, k)))*iareaT(i, j)
         end do
      end if
   end subroutine tracer_hdiff_one_impl