compute_fv_mom6_reconstruct_impl Subroutine

   pure subroutine compute_fv_mom6_reconstruct_impl(h_layer, hS, hT, b, eos, &
                                                    S_t, S_b, T_t, T_b, &
                                                    conc_T, conc_S, &
                                                    e_face, pa, intz_dpa, &
                                                    intx_pa, inty_pa, &
                                                    intx_dpa, inty_dpa, &
                                                    dpdx_face, dpdy_face, &
                                                    rho0, rho_ref, h_neglect, &
                                                    gfs_scale, recon_scheme, &
                                                    p_top, p_top_in_bc, &
                                                    idxCu, idyCv, nx, ny, nz)
      !! FV_MOM6 pressure-gradient with in-layer T/S reconstruction.
      !!
      !! Same Pass 3-5 face assembly as `compute_fv_mom6_impl`, but the
      !! two integrals that assembly consumes are BOTH taken from the
      !! reconstructed sub-layer T/S profile rather than a layer mean:
      !!
      !!   * Pass 1 replaces the PCM `dpa(k)` / `intz_dpa(k)` with the
      !!     5-point VERTICAL Boole quadrature of the monotone PLM/PPM
      !!     profile (edges from Pass 0) — the side integrals of the
      !!     control volume.
      !!   * Pass 2 (per face) replaces the two-column trapezoid
      !!     `0.5*(dpa_L + dpa_R)` with the 5-point HORIZONTAL Boole
      !!     quadrature `boole_dpa_face` — the top/bottom (tilted) edges.
      !!
      !! Passes 0-2 each run one GPU thread per CELL (3-D `do concurrent`
      !! over k, j, i, plus a cheap per-column scan for the `pa` / `intx_pa`
      !! / `inty_pa` recurrences) with every per-EOS helper inlined: the
      !! same operations in the same order as the column-serial form they
      !! replaced, so bit-identical to it.  Global 1-degree PPM, 5 days, one
      !! V100: `ocean_pgf` 5.97 -> 3.45 s under Wright, 12.42 -> 6.55 s under
      !! Roquet (`[stats]` identical to the digit).
      !!
      !! Both are required for the defining property: with a linear EOS
      !! and T/S linear in z, the PGF then vanishes to round-off for ANY
      !! layer geometry (Adcroft, Hallberg & Harrison 2008; Yung,
      !! Hallberg, Adcroft & Morrison 2026 §2.4).  Correcting the vertical
      !! integral alone leaves the horizontal trapezoid's curvature
      !! residual `g*(-drho/dz)*Delta_e^2/12` at every tilted interface,
      !! which is the sigma "second-kind" pressure-gradient error.
      !!
      !! `mass_weight` (hWght blend) is NOT applied here: it needs a
      !! per-cell density, whereas reconstruction works on column T/S
      !! edges.  Boundary layers take the linear-exact one-sided edge pair
      !! in the edge helper (`boundary_edges_linear`).
      !!
      !! `p_top_in_bc` injects the top load into the SAME Pass-1 surface
      !! BC as the PCM twin, and the Theorem in `compute_fv_mom6_impl`
      !! carries over verbatim: the reconstruction only changes `dpa` /
      !! `intz_dpa`, never the `pa(nz+1)` seed or the `intx_pa`
      !! recurrence, so a depth-uniform `p_top` still perturbs every
      !! layer's `PFu` by the same `−(1/ρ₀)∇p_top`. NOTE that this branch
      !! builds its OWN in-layer EOS pressure inside
      !! `boole_dpa_intz_layer` (`p = −g·ρ₀·z` from the surface-relative
      !! interface height) and that one is NOT offset by `p_top` — which
      !! is exactly why `validate_config` refuses
      !! `&ocean_psurf_nml in_eos` together with
      !! `reconstruct_for_pressure`. The BC injection here is a PRESSURE
      !! boundary condition, not an EOS argument; the two are independent
      !! seams.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in)    :: h_layer(nx, ny, nz)
      real(wp), intent(in)    :: hS(nx, ny, nz)
         !! Salinity * thickness (PSU*m) — layer-mean S = hS / h.
      real(wp), intent(in)    :: hT(nx, ny, nz)
         !! Temperature * thickness (degC*m) — layer-mean T = hT / h.
      real(wp), intent(in)    :: b(nx, ny)
      type(eos_t), intent(in) :: eos
      real(wp), intent(inout) :: S_t(nx, ny, nz), S_b(nx, ny, nz)
      real(wp), intent(inout) :: T_t(nx, ny, nz), T_b(nx, ny, nz)
      real(wp), intent(inout) :: conc_T(nx, ny, nz), conc_S(nx, ny, nz)
         !! Layer-mean T / S as every pass below reads them (Pass C).
      real(wp), intent(inout) :: e_face(nx, ny, nz + 1)
      real(wp), intent(inout) :: pa(nx, ny, nz + 1)
      real(wp), intent(inout) :: intz_dpa(nx, ny, nz)
      real(wp), intent(inout) :: intx_pa(nx + 1, ny, nz + 1)
      real(wp), intent(inout) :: inty_pa(nx, ny + 1, nz + 1)
      real(wp), intent(inout) :: intx_dpa(nx + 1, ny, nz)
      real(wp), intent(inout) :: inty_dpa(nx, ny + 1, nz)
      real(wp), intent(inout) :: dpdx_face(nx + 1, ny, nz)
      real(wp), intent(inout) :: dpdy_face(nx, ny + 1, nz)
      real(wp), intent(in)    :: rho0, rho_ref, h_neglect, gfs_scale
      integer, intent(in)    :: recon_scheme
      real(wp), intent(in)    :: p_top(nx, ny)
         !! Top-of-column pressure (Pa, `>= 0`), `multilayer_state_t%p_top`.
      logical, intent(in)    :: p_top_in_bc
         !! Add `p_top` to the Pass-1 surface BC (`.false.` ⇒ bit-identical).
      real(wp), intent(in)    :: idxCu(nx + 1, ny), idyCv(nx, ny + 1)

      integer  :: i, j, k
      real(wp) :: inv_rho0, dpa_kk, intz_kk
      real(wp) :: h_L, h_R, e_bot_L, e_bot_R
      real(wp) :: t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R
      real(wp) :: pa_h_intz_L, pa_h_intz_R, numer, denom
      real(wp) :: dM_coeff, ddM_dx, ddM_dy
      logical  :: parabolic
      integer  :: km2, km1, kp1, kp2, k_top_c, k_don_c
      integer  :: eos_variant

      inv_rho0 = 1.0_wp/rho0
      parabolic = (recon_scheme == PGF_RECON_PPM)
      eos_variant = eos%variant

      ! ---- Pass C (per column): the layer-mean T/S every pass reads ----
      ! `hT/h` on a live layer and, on a vanished one, its I1′ DONOR's:
      ! the nearest live layer above it, or for a run of fillers reaching
      ! the top of the column the topmost live layer; 0 in a column with
      ! no live layer.  The per-column form of `rdb_vl_column_conc`, read
      ! off the donor so no near-zero thickness is ever a divisor.  MOM6
      ! carries T/S as concentrations, so its vanished layers hold the
      ! remapped value its `int_density_dz_*` reads; `c_live` is that value
      ! here.  NOT the floored `hT/max(h, H_VANISHED)` this replaced: that
      ! is `h/H_VANISHED` of the truth on a filler (2/3 at the default
      ! `zstar_h_min = 1e-4 m`), harmless in the vertical `pa` stack where
      ! it multiplies the filler's own thickness, but the cross-face Boole
      ! integral interpolates T/S over the INTERPOLATED -- live --
      ! thickness and the PLM/PPM stencil reads its neighbours, so at an
      ! OPEN z-like step (`zstar`, closed-faces-off `z_fixed`, a `z_fixed`
      ! cell whose liveness flipped with eta under a static closed-face
      ! mask) it integrated the wrong salinity over tens of metres of live
      ! water: 2.9e-3 m/s^2 at rest on a live|filler face against 1.4e-6
      ! (`test_ocean_pgf_insitu :: open_step_filler_faces_*`).  One O(nz)
      ! sweep per column, the shape of Pass 1a: a per-cell donor walk made
      ! `ocean_pgf` 4x slower on the global 1-degree grid (long bed-filler
      ! runs).  A live layer reads `hT/h` exactly as before (bit-identical
      ! on a column without fillers).
      do concurrent(j=1:ny, i=1:nx) local(k, k_top_c, k_don_c)
         k_top_c = 0
         do k = nz, 1, -1
            if (rdb_vl_is_live(h_layer(i, j, k))) then
               k_top_c = k
               exit
            end if
         end do
         if (k_top_c == 0) then
            do k = 1, nz
               conc_T(i, j, k) = 0.0_wp
               conc_S(i, j, k) = 0.0_wp
            end do
         else
            k_don_c = k_top_c
            do k = nz, 1, -1
               if (rdb_vl_is_live(h_layer(i, j, k))) k_don_c = k
               conc_T(i, j, k) = rdb_vl_conc(hT(i, j, k_don_c), h_layer(i, j, k_don_c))
               conc_S(i, j, k) = rdb_vl_conc(hS(i, j, k_don_c), h_layer(i, j, k_don_c))
            end do
         end if
      end do

      ! ---- Pass 0: PLM/PPM T/S edge values, one thread per cell ----
      ! A layer's edges come from its own short vertical stencil of layer
      ! means (k-1..k+1 for PLM, k-2..k+2 for PPM), so this is a 3-D
      ! `do concurrent`; the per-column form built seven NZ_STACK_MAX
      ! stacks per thread in device local memory.  The means are Pass C's
      ! `conc_T`/`conc_S` (the I1′ donor's on a vanished layer).  Stencil
      ! indices outside the column are clamped into it; the boundary
      ! branches that would read them do not.  Done as its own pass so the
      ! quadrature passes below read clean edge arrays.
      do concurrent(k=1:nz, j=1:ny, i=1:nx) local(km2, km1, kp1, kp2)
         km2 = max(k - 2, 1)
         km1 = max(k - 1, 1)
         kp1 = min(k + 1, nz)
         kp2 = min(k + 2, nz)
         if (parabolic) then
            call ppm_edges_layer(k, nz, h_layer(i, j, km2), h_layer(i, j, km1), &
                                 h_layer(i, j, k), h_layer(i, j, kp1), h_layer(i, j, kp2), &
                                 conc_S(i, j, km2), &
                                 conc_S(i, j, km1), &
                                 conc_S(i, j, k), &
                                 conc_S(i, j, kp1), &
                                 conc_S(i, j, kp2), &
                                 conc_T(i, j, km2), &
                                 conc_T(i, j, km1), &
                                 conc_T(i, j, k), &
                                 conc_T(i, j, kp1), &
                                 conc_T(i, j, kp2), &
                                 S_t(i, j, k), S_b(i, j, k), T_t(i, j, k), T_b(i, j, k))
         else
            call plm_edges_layer(k, nz, h_layer(i, j, km1), h_layer(i, j, k), h_layer(i, j, kp1), &
                                 conc_S(i, j, km1), &
                                 conc_S(i, j, k), &
                                 conc_S(i, j, kp1), &
                                 S_t(i, j, k), S_b(i, j, k))
            call plm_edges_layer(k, nz, h_layer(i, j, km1), h_layer(i, j, k), h_layer(i, j, kp1), &
                                 conc_T(i, j, km1), &
                                 conc_T(i, j, k), &
                                 conc_T(i, j, kp1), &
                                 T_t(i, j, k), T_b(i, j, k))
         end if
      end do

      ! ---- Pass 1: e_face, pa, intz_dpa via Boole quadrature ----
      ! e_top = e_face(k+1) is the shallower interface of layer k.  The
      ! reconstructed dpa(k) marches the pa stack; intz_dpa(k) is the
      ! first-moment piece.  Both replace the PCM forms.
      !
      ! Pass 1a (per column): interface heights and the surface seed of the
      ! pressure-anomaly stack.  Pass 1b (3-D `do concurrent` over k, j, i):
      ! every layer's `dpa` / `intz_dpa`, `dpa` parked in `pa(i, j, k)`.
      ! Pass 1c (per column): the stack sum, top down, in place.  The same
      ! operations in the same order as the column-serial march, so
      ! bit-identical to it, but one GPU thread per CELL -- the per-column
      ! form ran one thread per column through nz layers x 5 EOS
      ! evaluations (115k threads on the global 1-degree grid).  Pass 2
      ! likewise: 3-D face integrals, then a per-column scan.
      !
      ! Passes 1b and 2 come in one loop copy per EOS, selected here, outside
      ! the loops (same sub-points, weights and summation order in every
      ! copy): Wright calls its Boole twins (density inline, no `eos_t`
      ! handle); Roquet calls this module's twins `roquet_recon_dpa_intz` /
      ! `roquet_recon_dpa_face`, whose SpV value comes from the included
      ! `rdb_roquet_spv.inc` and is inlined into the kernel (the generic
      ! chain's out-of-line `eos_density_point` calls were this path's
      ! device cost); anything else (the linear EOS) takes the generic
      ! `eos_density_point` rule.
      do concurrent(j=1:ny, i=1:nx) local(k)
         e_face(i, j, 1) = -b(i, j)
         do k = 1, nz
            e_face(i, j, k + 1) = e_face(i, j, k) + h_layer(i, j, k)
         end do
         if (p_top_in_bc) then
            pa(i, j, nz + 1) = rho_ref*GRAVITY*e_face(i, j, nz + 1) + p_top(i, j)
         else
            pa(i, j, nz + 1) = rho_ref*GRAVITY*e_face(i, j, nz + 1)
         end if
      end do
      if (eos_variant == EOS_VARIANT_WRIGHT_97) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx) local(dpa_kk, intz_kk)
            call boole_dpa_intz_layer_wright(rho0, rho_ref, &
                                             e_face(i, j, k + 1), h_layer(i, j, k), &
                                             T_t(i, j, k), T_b(i, j, k), &
                                             conc_T(i, j, k), &
                                             S_t(i, j, k), S_b(i, j, k), &
                                             conc_S(i, j, k), &
                                             parabolic, dpa_kk, intz_kk)
            pa(i, j, k) = dpa_kk
            intz_dpa(i, j, k) = intz_kk
         end do
      else if (eos_variant == EOS_VARIANT_ROQUET_SPV) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx) local(dpa_kk, intz_kk)
            call roquet_recon_dpa_intz(rho0, rho_ref, &
                                       e_face(i, j, k + 1), h_layer(i, j, k), &
                                       T_t(i, j, k), T_b(i, j, k), &
                                       conc_T(i, j, k), &
                                       S_t(i, j, k), S_b(i, j, k), &
                                       conc_S(i, j, k), &
                                       parabolic, dpa_kk, intz_kk)
            pa(i, j, k) = dpa_kk
            intz_dpa(i, j, k) = intz_kk
         end do
      else
         do concurrent(k=1:nz, j=1:ny, i=1:nx) local(dpa_kk, intz_kk)
            call boole_dpa_intz_layer(eos, rho0, rho_ref, &
                                      e_face(i, j, k + 1), h_layer(i, j, k), &
                                      T_t(i, j, k), T_b(i, j, k), &
                                      conc_T(i, j, k), &
                                      S_t(i, j, k), S_b(i, j, k), &
                                      conc_S(i, j, k), &
                                      parabolic, dpa_kk, intz_kk)
            pa(i, j, k) = dpa_kk
            intz_dpa(i, j, k) = intz_kk
         end do
      end if
      do concurrent(j=1:ny, i=1:nx) local(k)
         do k = nz, 1, -1
            pa(i, j, k) = pa(i, j, k + 1) + pa(i, j, k)
         end do
      end do

      ! ---- Pass 2a: u-face horizontal integrals ----
      ! The along-face mean of the layer pressure increment, by the 5-point
      ! cross-face Boole quadrature of `boole_dpa_face` (sub-columns at the
      ! INTERPOLATED interface height with interpolated T/S).  The
      ! two-column trapezoid `0.5*(dpa_L + dpa_R)` this replaces is exact
      ! only for a pressure linear in x along the edge; under a tilted
      ! interface it leaves the sigma second-kind curvature residual at
      ! every interface — see the `boole_dpa_face` docstring.
      if (eos_variant == EOS_VARIANT_WRIGHT_97) then
         do concurrent(k=1:nz, j=1:ny, i=2:nx) &
            local(dpa_kk, t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R)
            t_m_L = conc_T(i - 1, j, k)
            t_m_R = conc_T(i, j, k)
            s_m_L = conc_S(i - 1, j, k)
            s_m_R = conc_S(i, j, k)
            dpa_L = pa(i - 1, j, k) - pa(i - 1, j, k + 1)
            dpa_R = pa(i, j, k) - pa(i, j, k + 1)
            call boole_dpa_face_wright(rho0, rho_ref, &
                                       e_face(i - 1, j, k + 1), e_face(i, j, k + 1), &
                                       h_layer(i - 1, j, k), h_layer(i, j, k), &
                                       T_t(i - 1, j, k), T_b(i - 1, j, k), t_m_L, &
                                       T_t(i, j, k), T_b(i, j, k), t_m_R, &
                                       S_t(i - 1, j, k), S_b(i - 1, j, k), s_m_L, &
                                       S_t(i, j, k), S_b(i, j, k), s_m_R, &
                                       dpa_L, dpa_R, parabolic, dpa_kk)
            intx_dpa(i, j, k) = dpa_kk
         end do
      else if (eos_variant == EOS_VARIANT_ROQUET_SPV) then
         do concurrent(k=1:nz, j=1:ny, i=2:nx) &
            local(dpa_kk, t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R)
            t_m_L = conc_T(i - 1, j, k)
            t_m_R = conc_T(i, j, k)
            s_m_L = conc_S(i - 1, j, k)
            s_m_R = conc_S(i, j, k)
            dpa_L = pa(i - 1, j, k) - pa(i - 1, j, k + 1)
            dpa_R = pa(i, j, k) - pa(i, j, k + 1)
            call roquet_recon_dpa_face(rho0, rho_ref, &
                                       e_face(i - 1, j, k + 1), e_face(i, j, k + 1), &
                                       h_layer(i - 1, j, k), h_layer(i, j, k), &
                                       T_t(i - 1, j, k), T_b(i - 1, j, k), t_m_L, &
                                       T_t(i, j, k), T_b(i, j, k), t_m_R, &
                                       S_t(i - 1, j, k), S_b(i - 1, j, k), s_m_L, &
                                       S_t(i, j, k), S_b(i, j, k), s_m_R, &
                                       dpa_L, dpa_R, parabolic, dpa_kk)
            intx_dpa(i, j, k) = dpa_kk
         end do
      else
         do concurrent(k=1:nz, j=1:ny, i=2:nx) &
            local(dpa_kk, t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R)
            t_m_L = conc_T(i - 1, j, k)
            t_m_R = conc_T(i, j, k)
            s_m_L = conc_S(i - 1, j, k)
            s_m_R = conc_S(i, j, k)
            dpa_L = pa(i - 1, j, k) - pa(i - 1, j, k + 1)
            dpa_R = pa(i, j, k) - pa(i, j, k + 1)
            call boole_dpa_face(eos, rho0, rho_ref, &
                                e_face(i - 1, j, k + 1), e_face(i, j, k + 1), &
                                h_layer(i - 1, j, k), h_layer(i, j, k), &
                                T_t(i - 1, j, k), T_b(i - 1, j, k), t_m_L, &
                                T_t(i, j, k), T_b(i, j, k), t_m_R, &
                                S_t(i - 1, j, k), S_b(i - 1, j, k), s_m_L, &
                                S_t(i, j, k), S_b(i, j, k), s_m_R, &
                                dpa_L, dpa_R, parabolic, dpa_kk)
            intx_dpa(i, j, k) = dpa_kk
         end do
      end if
      ! Column scan of the face integrals (cheap; the EOS work above is
      ! 3-D parallel).
      do concurrent(j=1:ny, i=2:nx) local(k)
         intx_pa(i, j, nz + 1) = 0.5_wp*(pa(i - 1, j, nz + 1) + pa(i, j, nz + 1))
         do k = nz, 1, -1
            intx_pa(i, j, k) = intx_pa(i, j, k + 1) + intx_dpa(i, j, k)
         end do
      end do
      do concurrent(k=1:nz, j=1:ny)
         intx_dpa(1, j, k) = 0.0_wp
         intx_dpa(nx + 1, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz + 1, j=1:ny)
         intx_pa(1, j, k) = 0.0_wp
         intx_pa(nx + 1, j, k) = 0.0_wp
      end do

      ! ---- Pass 2b: v-face horizontal integrals ----
      if (eos_variant == EOS_VARIANT_WRIGHT_97) then
         do concurrent(k=1:nz, j=2:ny, i=1:nx) &
            local(dpa_kk, t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R)
            t_m_L = conc_T(i, j - 1, k)
            t_m_R = conc_T(i, j, k)
            s_m_L = conc_S(i, j - 1, k)
            s_m_R = conc_S(i, j, k)
            dpa_L = pa(i, j - 1, k) - pa(i, j - 1, k + 1)
            dpa_R = pa(i, j, k) - pa(i, j, k + 1)
            call boole_dpa_face_wright(rho0, rho_ref, &
                                       e_face(i, j - 1, k + 1), e_face(i, j, k + 1), &
                                       h_layer(i, j - 1, k), h_layer(i, j, k), &
                                       T_t(i, j - 1, k), T_b(i, j - 1, k), t_m_L, &
                                       T_t(i, j, k), T_b(i, j, k), t_m_R, &
                                       S_t(i, j - 1, k), S_b(i, j - 1, k), s_m_L, &
                                       S_t(i, j, k), S_b(i, j, k), s_m_R, &
                                       dpa_L, dpa_R, parabolic, dpa_kk)
            inty_dpa(i, j, k) = dpa_kk
         end do
      else if (eos_variant == EOS_VARIANT_ROQUET_SPV) then
         do concurrent(k=1:nz, j=2:ny, i=1:nx) &
            local(dpa_kk, t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R)
            t_m_L = conc_T(i, j - 1, k)
            t_m_R = conc_T(i, j, k)
            s_m_L = conc_S(i, j - 1, k)
            s_m_R = conc_S(i, j, k)
            dpa_L = pa(i, j - 1, k) - pa(i, j - 1, k + 1)
            dpa_R = pa(i, j, k) - pa(i, j, k + 1)
            call roquet_recon_dpa_face(rho0, rho_ref, &
                                       e_face(i, j - 1, k + 1), e_face(i, j, k + 1), &
                                       h_layer(i, j - 1, k), h_layer(i, j, k), &
                                       T_t(i, j - 1, k), T_b(i, j - 1, k), t_m_L, &
                                       T_t(i, j, k), T_b(i, j, k), t_m_R, &
                                       S_t(i, j - 1, k), S_b(i, j - 1, k), s_m_L, &
                                       S_t(i, j, k), S_b(i, j, k), s_m_R, &
                                       dpa_L, dpa_R, parabolic, dpa_kk)
            inty_dpa(i, j, k) = dpa_kk
         end do
      else
         do concurrent(k=1:nz, j=2:ny, i=1:nx) &
            local(dpa_kk, t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R)
            t_m_L = conc_T(i, j - 1, k)
            t_m_R = conc_T(i, j, k)
            s_m_L = conc_S(i, j - 1, k)
            s_m_R = conc_S(i, j, k)
            dpa_L = pa(i, j - 1, k) - pa(i, j - 1, k + 1)
            dpa_R = pa(i, j, k) - pa(i, j, k + 1)
            call boole_dpa_face(eos, rho0, rho_ref, &
                                e_face(i, j - 1, k + 1), e_face(i, j, k + 1), &
                                h_layer(i, j - 1, k), h_layer(i, j, k), &
                                T_t(i, j - 1, k), T_b(i, j - 1, k), t_m_L, &
                                T_t(i, j, k), T_b(i, j, k), t_m_R, &
                                S_t(i, j - 1, k), S_b(i, j - 1, k), s_m_L, &
                                S_t(i, j, k), S_b(i, j, k), s_m_R, &
                                dpa_L, dpa_R, parabolic, dpa_kk)
            inty_dpa(i, j, k) = dpa_kk
         end do
      end if
      do concurrent(j=2:ny, i=1:nx) local(k)
         inty_pa(i, j, nz + 1) = 0.5_wp*(pa(i, j - 1, nz + 1) + pa(i, j, nz + 1))
         do k = nz, 1, -1
            inty_pa(i, j, k) = inty_pa(i, j, k + 1) + inty_dpa(i, j, k)
         end do
      end do
      do concurrent(k=1:nz, i=1:nx)
         inty_dpa(i, 1, k) = 0.0_wp
         inty_dpa(i, ny + 1, k) = 0.0_wp
      end do
      do concurrent(k=1:nz + 1, i=1:nx)
         inty_pa(i, 1, k) = 0.0_wp
         inty_pa(i, ny + 1, k) = 0.0_wp
      end do

      ! ---- Pass 3: PFu assembly (identical to compute_fv_mom6_impl) ----
      do concurrent(k=1:nz, j=1:ny, i=2:nx) &
         local(h_L, h_R, e_bot_L, e_bot_R, pa_h_intz_L, pa_h_intz_R, numer, denom)
         h_L = h_layer(i - 1, j, k)
         h_R = h_layer(i, j, k)
         e_bot_L = e_face(i - 1, j, k)
         e_bot_R = e_face(i, j, k)
         pa_h_intz_L = pa(i - 1, j, k + 1)*h_L + intz_dpa(i - 1, j, k)
         pa_h_intz_R = pa(i, j, k + 1)*h_R + intz_dpa(i, j, k)
         numer = (pa_h_intz_L - pa_h_intz_R) &
                 + (h_R - h_L)*intx_pa(i, j, k + 1) &
                 - (e_bot_R - e_bot_L)*intx_dpa(i, j, k)
         denom = h_L + h_R + h_neglect
         dpdx_face(i, j, k) = numer*(2.0_wp*inv_rho0*idxCu(i, j))/denom
      end do
      do concurrent(k=1:nz, j=1:ny)
         dpdx_face(1, j, k) = 0.0_wp
         dpdx_face(nx + 1, j, k) = 0.0_wp
      end do

      ! ---- Pass 4: PFv assembly ----
      do concurrent(k=1:nz, j=2:ny, i=1:nx) &
         local(h_L, h_R, e_bot_L, e_bot_R, pa_h_intz_L, pa_h_intz_R, numer, denom)
         h_L = h_layer(i, j - 1, k)
         h_R = h_layer(i, j, k)
         e_bot_L = e_face(i, j - 1, k)
         e_bot_R = e_face(i, j, k)
         pa_h_intz_L = pa(i, j - 1, k + 1)*h_L + intz_dpa(i, j - 1, k)
         pa_h_intz_R = pa(i, j, k + 1)*h_R + intz_dpa(i, j, k)
         numer = (pa_h_intz_L - pa_h_intz_R) &
                 + (h_R - h_L)*inty_pa(i, j, k + 1) &
                 - (e_bot_R - e_bot_L)*inty_dpa(i, j, k)
         denom = h_L + h_R + h_neglect
         dpdy_face(i, j, k) = numer*(2.0_wp*inv_rho0*idyCv(i, j))/denom
      end do
      do concurrent(k=1:nz, i=1:nx)
         dpdy_face(i, 1, k) = 0.0_wp
         dpdy_face(i, ny + 1, k) = 0.0_wp
      end do

      ! ---- Pass 5: Montgomery dM correction (MOM6 GFS_scale) ----
      ! rho_surf for the dM term is the reconstructed top-edge density of
      ! the surface layer; here we reuse the layer-mean surface density
      ! recovered from dpa(nz) = pa(nz) - pa(nz+1) divided by g*h, which
      ! equals (rho_surf - rho_ref).  Keep the same depth-independent form.
      if (gfs_scale < 1.0_wp - 1.0e-12_wp) then
         dM_coeff = (gfs_scale - 1.0_wp)*GRAVITY*inv_rho0
         do concurrent(k=1:nz, j=1:ny, i=2:nx) local(ddM_dx)
            ddM_dx = dM_coeff*(recon_rho_surf(pa(i, j, nz), pa(i, j, nz + 1), &
                                              h_layer(i, j, nz), rho_ref) &
                               *e_face(i, j, nz + 1) &
                               - recon_rho_surf(pa(i - 1, j, nz), pa(i - 1, j, nz + 1), &
                                                h_layer(i - 1, j, nz), rho_ref) &
                               *e_face(i - 1, j, nz + 1))*idxCu(i, j)
            dpdx_face(i, j, k) = dpdx_face(i, j, k) - ddM_dx
         end do
         do concurrent(k=1:nz, j=2:ny, i=1:nx) local(ddM_dy)
            ddM_dy = dM_coeff*(recon_rho_surf(pa(i, j, nz), pa(i, j, nz + 1), &
                                              h_layer(i, j, nz), rho_ref) &
                               *e_face(i, j, nz + 1) &
                               - recon_rho_surf(pa(i, j - 1, nz), pa(i, j - 1, nz + 1), &
                                                h_layer(i, j - 1, nz), rho_ref) &
                               *e_face(i, j - 1, nz + 1))*idyCv(i, j)
            dpdy_face(i, j, k) = dpdy_face(i, j, k) - ddM_dy
         end do
      end if
   end subroutine compute_fv_mom6_reconstruct_impl