compute_fv_mom6_insitu_pcm_impl Subroutine

   pure subroutine compute_fv_mom6_insitu_pcm_impl(h_layer, hS, hT, 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, mass_weight, &
                                                   eos_variant, &
                                                   p_top, p_top_in_bc, &
                                                   idxCu, idyCv, nx, ny, nz)
      !! FV_MOM6 pressure gradient, constant-by-layer (PCM) T/S, density at
      !! the IN-SITU pressure — MOM6 `PressureForce_FV_Bouss` with
      !! `RECONSTRUCT_FOR_PRESSURE = False` (`int_density_dz_generic_pcm`).
      !!
      !! The PCM twin `compute_fv_mom6_impl` integrates `ms%rho_layer`, a
      !! POTENTIAL density at the one horizontally uniform `p_ref`.  Its
      !! horizontal difference at depth is then the difference at the
      !! REFERENCE pressure, not at the local one: the thermal expansion
      !! coefficient roughly doubles between the surface and 4000 dbar
      !! (thermobaricity), so with `p_ref = 0` the deep baroclinic
      !! pressure gradient — the bottom-pressure gradient that forces the
      !! barotropic mode over topography — is systematically too weak.
      !! On the global 1-degree WOA13 spin-up it held Drake Passage at
      !! ~80 Sv where MOM6 on the same protocol adjusts to ~155 Sv, and
      !! the transport tracked `p_ref` (0 / 2000 / 4000 dbar: 83 / 143 /
      !! 203 Sv) — the tell of a reference-pressure artefact.
      !!
      !! Here every density is `EOS(T, S, p = −g·rho0·z)` at the point it
      !! is used, integrated (Roquet; Wright takes the closed form below)
      !! by the same 5-point Boole rules as the reconstruction branch, with
      !! the sub-layer profile flat:
      !!
      !!   * Pass 1 (per column): `dpa(k)`, `intz_dpa(k)` from
      !!     `boole_dpa_intz_layer` with top = bottom = mean T/S.
      !!   * Pass 2 (per face): `intx_dpa` / `inty_dpa` from
      !!     `boole_dpa_face_pcm` — end points are the columns' own `dpa`,
      !!     the three interior sub-columns interpolate `z` linearly and
      !!     T/S with MOM6's near-bottom mass weighting (`hWght`, the same
      !!     measure and blend as `compute_fv_mom6_impl`) when
      !!     `mass_weight`.
      !!   * Passes 3–5: the face assembly, identical to the other two
      !!     FV_MOM6 branches.
      !!
      !! Under Wright (`wright_analytic`) Passes 1–2 replace each vertical
      !! Boole rule by the closed-form integral `wright_pcm_dpa_intz` (MOM6
      !! `int_density_dz_wright`), keeping the 5-point cross-face Boole
      !! rule and the same sub-columns (`wright_pcm_dpa_face`).  The
      !! vertical Boole rule it replaces is accurate to
      !! `(g·rho0·dz/(p + p0 + lambda/alpha0))^6` — round-off for any
      !! realistic layer — so the answers move at round-off, but each layer
      !! costs one polynomial evaluation instead of 5 generic-EOS calls and
      !! each face 3 instead of 15.  Measured on the global 1° run (5
      !! days, one V100): `ocean_pgf` 10.26 s → 1.90 s, against 1.09 s for
      !! `insitu_density = .false.`.
      !!
      !! Under Roquet the vertical rule stays 5-point Boole (no closed form
      !! for `int dz/SV(p)`), factored: `roquet_pcm_dpa_intz` evaluates the
      !! (T, S) part of the SpV polynomial once per sub-column and only the
      !! pressure Horner per point — the same five densities as
      !! `boole_dpa_intz_layer`, so answers are unchanged (to the digit on
      !! the global 1° run's `[stats]` and En).  `ocean_pgf` 10.73 s →
      !! 3.64 s on the same run, → 2.55 s (1.36x Wright) with the SpV value
      !! from the module-local include `rdb_roquet_spv.inc`: as a call into
      !! `rdb_eos` the (T, S) part was NOT inlined on the device (a real
      !! `call`, its four results through the stack, 130+ registers in the
      !! face kernels against 88 inlined).
      !!
      !! Pass 1 and Pass 2 each run their integrals as a 3-D `do concurrent`
      !! over (k, j, i) — every layer's and every face's integral is
      !! independent — followed by a cheap per-column scan for the `pa` /
      !! `intx_pa` / `inty_pa` recurrences (Pass 1 parks `dpa` in `pa(k)`
      !! and sums it in place).  Same operations in the same order, so
      !! bit-identical to the column-serial form; one GPU thread per CELL
      !! instead of per column.
      !!
      !! The trapezoid `0.5·(dpa_L + dpa_R)` of the potential-density twin
      !! is NOT kept: an in-situ density carries the compressibility
      !! gradient (`~4.4e-3 kg m⁻⁴`), and the trapezoid's curvature
      !! residual `g·(∂ρ/∂z)·Δe²/12` at a tilted interface (a partial-cell
      !! bed step) would be of the size of the signal.
      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)
      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
      logical, intent(in)    :: mass_weight
         !! MOM6 `MASS_WEIGHT_IN_PRESSURE_GRADIENT` (near-bottom `hWght`).
      integer, intent(in)    :: eos_variant
         !! `eos%variant`, selected ONCE here, outside the loops -- each
         !! variant has its own loop copies with its integrals inlined, no
         !! generic `eos_t` dispatch in the hot loop.  `use_insitu_pcm`
         !! admits exactly two:
         !!   * `EOS_VARIANT_WRIGHT_97`: the ANALYTIC Wright layer integral
         !!     (`wright_pcm_dpa_intz` / `wright_pcm_dpa_face`, MOM6
         !!     `int_density_dz_wright`).
         !!   * `EOS_VARIANT_ROQUET_SPV` (the `else` copies): the 5-point
         !!     Boole rule of MOM6 `int_density_dz_generic_pcm`, with the
         !!     (T, S) part of the EOS evaluated once per sub-column
         !!     (`roquet_pcm_dpa_intz` / `roquet_pcm_dpa_face`).
      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.` ⇒ the plain
         !! `rho_ref·g·eta` seed).
      real(wp), intent(in)    :: idxCu(nx + 1, ny), idyCv(nx, ny + 1)

      integer  :: i, j, k
      real(wp) :: inv_rho0, dpa_kk, intz_kk, t_m, s_m
      real(wp) :: t_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R
      real(wp) :: hwt_ll, hwt_lr, hwt_rr, hwt_rl
      real(wp) :: h_L, h_R, e_bot_L, e_bot_R
      real(wp) :: pa_h_intz_L, pa_h_intz_R, numer, denom
      real(wp) :: dM_coeff, ddM_dx, ddM_dy
      logical  :: wright_analytic
      integer  :: k_top_c, k_don_c

      inv_rho0 = 1.0_wp/rho0
      wright_analytic = (eos_variant == EOS_VARIANT_WRIGHT_97)

      ! ---- 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 1: per-column e_face, pa, intz_dpa (in-situ) ----
      ! Pass 1a (per column): interface heights and the surface seed of
      ! the pressure-anomaly stack.
      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
      ! Pass 1b (3-D parallel): every layer's `dpa` / `intz_dpa`.  `dpa`
      ! is parked in `pa(i, j, k)` and summed by Pass 1c.  Layer by layer
      ! the EOS work is independent; only the stack sum is a recurrence, so
      ! the expensive part no longer runs one thread per COLUMN (115k
      ! threads, latency-bound on the global 1-degree grid) but one per
      ! CELL.
      ! Two loop copies (not a branch inside one kernel) so the analytic
      ! Wright kernel carries none of the Boole path's register pressure.
      if (wright_analytic) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx) local(dpa_kk, intz_kk, t_m, s_m)
            t_m = conc_T(i, j, k)
            s_m = conc_S(i, j, k)
            call wright_pcm_dpa_intz(t_m, s_m, e_face(i, j, k + 1), h_layer(i, j, k), &
                                     rho0, rho_ref, 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, t_m, s_m)
            t_m = conc_T(i, j, k)
            s_m = conc_S(i, j, k)
            call roquet_pcm_dpa_intz(t_m, s_m, e_face(i, j, k + 1), h_layer(i, j, k), &
                                     rho0, rho_ref, dpa_kk, intz_kk)
            pa(i, j, k) = dpa_kk
            intz_dpa(i, j, k) = intz_kk
         end do
      end if
      ! Pass 1c (per column): the pressure-anomaly stack, top down.
      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 ----
      if (wright_analytic) 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, hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            call fv_mom6_mass_weights(mass_weight, e_face(i - 1, j, 1), e_face(i, j, 1), &
                                      e_face(i - 1, j, k + 1), e_face(i, j, k + 1), &
                                      h_layer(i - 1, j, k), h_layer(i, j, k), h_neglect, &
                                      hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            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 wright_pcm_dpa_face(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_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R, &
                                     hwt_ll, hwt_lr, hwt_rr, hwt_rl, rho0, rho_ref, 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, hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            call fv_mom6_mass_weights(mass_weight, e_face(i - 1, j, 1), e_face(i, j, 1), &
                                      e_face(i - 1, j, k + 1), e_face(i, j, k + 1), &
                                      h_layer(i - 1, j, k), h_layer(i, j, k), h_neglect, &
                                      hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            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_pcm_dpa_face(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_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R, &
                                     hwt_ll, hwt_lr, hwt_rr, hwt_rl, rho0, rho_ref, 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 (wright_analytic) 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, hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            call fv_mom6_mass_weights(mass_weight, e_face(i, j - 1, 1), e_face(i, j, 1), &
                                      e_face(i, j - 1, k + 1), e_face(i, j, k + 1), &
                                      h_layer(i, j - 1, k), h_layer(i, j, k), h_neglect, &
                                      hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            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 wright_pcm_dpa_face(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_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R, &
                                     hwt_ll, hwt_lr, hwt_rr, hwt_rl, rho0, rho_ref, 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, hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            call fv_mom6_mass_weights(mass_weight, e_face(i, j - 1, 1), e_face(i, j, 1), &
                                      e_face(i, j - 1, k + 1), e_face(i, j, k + 1), &
                                      h_layer(i, j - 1, k), h_layer(i, j, k), h_neglect, &
                                      hwt_ll, hwt_lr, hwt_rr, hwt_rl)
            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_pcm_dpa_face(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_m_L, t_m_R, s_m_L, s_m_R, dpa_L, dpa_R, &
                                     hwt_ll, hwt_lr, hwt_rr, hwt_rl, rho0, rho_ref, dpa_kk)
            inty_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=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) ----
      ! Same depth-independent form as the reconstruction branch, with the
      ! surface layer's mean in-situ density recovered from its `dpa`.
      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_insitu_pcm_impl