ocean_pressure_force_compute Subroutine

   pure subroutine ocean_pressure_force_compute(grid, metrics, pgf, ms, eos)
      !! Compute the hydrostatic pressure-gradient acceleration at every
      !! C-grid face. Variants (see the OPGF_VARIANT_* / GPRIME / FV_MOM6
      !! constants for the per-variant formulas): MONT (layer-mean ρgh),
      !! FV_LITE (+ z-correction), FV_WRIGHT (+ in-situ Wright density),
      !! GPRIME, FV_MOM6.
      !!
      !! Pass layout: (1b) z_centre from h_layer (MONT + FV variants);
      !! then either (M1) the Montgomery column recursion and (M2/M3) its
      !! face passes, or (1) the column pressure sweep filling p_edge
      !! (+ rho_insitu for FV_WRIGHT) and (2/3) the FV east/north-face
      !! acceleration; finally (4) the grounded-layer gate.
      !!
      !! `ms%rho_layer` must be up to date — call the EOS kernel first.
      !! FV_WRIGHT needs `EOS_VARIANT_WRIGHT_97` for a good Picard seed.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
         !! Curvilinear horizontal metrics. The u-face gradient divides by
         !! `idxCu(i,j)`, the v-face by `idyCv(i,j)`. On uniform Cartesian
         !! `idxCu == 1/dx` bitwise (byte-identical to scalar inv_dx/inv_dy).
      type(ocean_pressure_force_t), intent(inout) :: pgf
      type(multilayer_state_t), intent(in) :: ms
      type(eos_t), intent(in), optional :: eos
         !! EOS handle — REQUIRED when `pgf%reconstruct_for_pressure` is
         !! on (the in-layer Boole quadrature evaluates the EOS at each
         !! sub-point).  Optional so the legacy PCM call sites (and the
         !! non-reconstruct variant tests) need not thread it through.

      integer :: i, j, k, nx, ny, nz
      real(wp) :: inv_rho0, g_over_rho0
      real(wp) :: p_centre_left, p_centre_right
      real(wp) :: p_centre_below, p_centre_above
      real(wp) :: rho_face, z_correction, z_running
      real(wp) :: h_l, h_r, e_l, e_r, z_eff, drho_star
      real(wp) :: vtol

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      inv_rho0 = 1.0_wp/pgf%rho0

      ! ---- gprime / reduced-gravity branch (Tier-1: NK = 2 only) ----
      if (pgf%variant == OPGF_VARIANT_GPRIME) then
         call compute_gprime_impl(ms%h_layer, pgf%b, &
                                  pgf%dpdx_face%data, pgf%dpdy_face%data, &
                                  pgf%gprime_gfs, pgf%gprime_gint, &
                                  metrics%idxCu, metrics%idyCv, nx, ny, nz)
         return
      end if

      ! ---- Pass 1b: per-column z_centre from h_layer ----
      ! z_centre(:, :, k) = physical z of the layer-k mid-depth,
      ! measured from the free surface downward (z=0 at surface,
      ! z negative below).  Computed by summing h_layer top-down so
      ! the surface is the reference point regardless of column
      ! bathymetry.  This is the load-bearing change that makes the
      ! FV_LITE / FV_WRIGHT Jacobian cancel the cross-bathymetry
      ! pressure gradient: a bed-relative z_centre would inject a
      ! spurious `-g·ρ·dH/dx` residual at every shelf-break face.
      !
      ! Hoisted AHEAD of the variant branches because it depends on nothing
      ! but `h_layer`.  FV_LITE / FV_WRIGHT consume it in Passes 2/3; MONT
      ! recovers the layer-TOP interface height `z_centre + h/2` from it for
      ! both the `M` recursion and the face `z_eff`.  FV_MOM6's own face
      ! assembly never touches it, so on that path it is filled ONLY to feed
      ! the Pass-4 grounded-layer gate — which is also exactly when
      ! `ocean_pressure_force_init`'s allocation gate provides the buffer.
      if (pgf%variant == OPGF_VARIANT_MONT .or. &
          pgf%variant == OPGF_VARIANT_FV_LITE .or. &
          pgf%variant == OPGF_VARIANT_FV_WRIGHT .or. &
          (pgf%variant == OPGF_VARIANT_FV_MOM6 .and. pgf%skip_nonoverlap)) then
         do concurrent(j=1:ny, i=1:nx) local(z_running)
            z_running = 0.0_wp
            do k = nz, 1, -1
               pgf%z_centre%data(i, j, k) = z_running - 0.5_wp*ms%h_layer(i, j, k)
               z_running = z_running - ms%h_layer(i, j, k)
            end do
         end do
      end if

      ! ---- MONT branch — Boussinesq Montgomery potential ----------------
      ! Pass M1 (per column): the vertical recursion that builds `M`.  Pass
      ! M2/M3 (per face): ONE horizontal difference of `M`, plus the
      ! horizontal-density term.  No pressure stack is built on this path.
      ! See the OPGF_VARIANT_MONT docstring for the derivation.
      if (pgf%variant == OPGF_VARIANT_MONT) then
         g_over_rho0 = GRAVITY*inv_rho0

         ! ---- Pass M1: Montgomery potential per column -------------------
         ! Seeded at the free surface: `e_edge(nz+1) = 0` (z_centre is
         ! surface-relative) and the surface pressure is zero, so
         ! `M(nz) = p/rho0 + rho_star(nz)*e_edge(nz+1)` is identically zero
         ! in EVERY column.  That is not a loss: the barotropic `-g*grad(eta)`
         ! it would otherwise carry is the BT substep's job, and adding it
         ! here would double-count it.  Accumulated straight into the array
         ! (no `local` reassigned in the k-loop), mirroring the Pass-1
         ! pressure sweep.
         do concurrent(j=1:ny, i=1:nx)
            pgf%mont_M%data(i, j, nz) = 0.0_wp
            do k = nz - 1, 1, -1
               ! `e_edge(k+1)` — the interface SHARED by layers k and k+1,
               ! i.e. the TOP of layer k — is where `p` and `z` agree between
               ! the two layers, so the whole jump in M is the jump in
               ! rho_star.  Recovered from the layer-k centre.
               pgf%mont_M%data(i, j, k) = pgf%mont_M%data(i, j, k + 1) + &
                                          g_over_rho0*(ms%rho_layer(i, j, k) - &
                                                       ms%rho_layer(i, j, k + 1))* &
                                          (pgf%z_centre%data(i, j, k) + &
                                           0.5_wp*ms%h_layer(i, j, k))
            end do
         end do

         ! ---- Pass M2: east-face acceleration ----------------------------
         ! `-dM/dx` plus the horizontal-density term `+ z_eff * d(rho_star)/dx`.
         ! `z_eff` is the thickness-weighted height at which `M`'s two
         ! column anchors are reconciled; on aligned columns (h_L == h_R) it
         ! is exactly the mean layer centre and the pair reduces
         ! ALGEBRAICALLY to FV_LITE.  `H_DIV_EPS` is pure 1/0 armour for a
         ! face between two fully-vanished layers (numerator is then zero
         ! too, so the face value is a clean zero).
         do concurrent(k=1:nz, j=1:ny, i=2:nx) &
            local(h_l, h_r, e_l, e_r, z_eff, drho_star)
            h_l = ms%h_layer(i - 1, j, k)
            h_r = ms%h_layer(i, j, k)
            e_l = pgf%z_centre%data(i - 1, j, k) + 0.5_wp*h_l
            e_r = pgf%z_centre%data(i, j, k) + 0.5_wp*h_r
            z_eff = (e_l*h_r + e_r*h_l - h_l*h_r)/(h_l + h_r + H_DIV_EPS)
            drho_star = g_over_rho0*(ms%rho_layer(i, j, k) - ms%rho_layer(i - 1, j, k))
            pgf%dpdx_face%data(i, j, k) = &
               (-(pgf%mont_M%data(i, j, k) - pgf%mont_M%data(i - 1, j, k)) &
                + drho_star*z_eff)*metrics%idxCu(i, j)
         end do
         do concurrent(k=1:nz, j=1:ny)
            pgf%dpdx_face%data(1, j, k) = 0.0_wp
            pgf%dpdx_face%data(nx + 1, j, k) = 0.0_wp
         end do

         ! ---- Pass M3: north-face acceleration ---------------------------
         do concurrent(k=1:nz, j=2:ny, i=1:nx) &
            local(h_l, h_r, e_l, e_r, z_eff, drho_star)
            h_l = ms%h_layer(i, j - 1, k)
            h_r = ms%h_layer(i, j, k)
            e_l = pgf%z_centre%data(i, j - 1, k) + 0.5_wp*h_l
            e_r = pgf%z_centre%data(i, j, k) + 0.5_wp*h_r
            z_eff = (e_l*h_r + e_r*h_l - h_l*h_r)/(h_l + h_r + H_DIV_EPS)
            drho_star = g_over_rho0*(ms%rho_layer(i, j, k) - ms%rho_layer(i, j - 1, k))
            pgf%dpdy_face%data(i, j, k) = &
               (-(pgf%mont_M%data(i, j, k) - pgf%mont_M%data(i, j - 1, k)) &
                + drho_star*z_eff)*metrics%idyCv(i, j)
         end do
         do concurrent(k=1:nz, i=1:nx)
            pgf%dpdy_face%data(i, 1, k) = 0.0_wp
            pgf%dpdy_face%data(i, ny + 1, k) = 0.0_wp
         end do

         ! ---- FV_MOM6 branch — faithful port of MOM6 PressureForce_FV_Bouss ----
         ! Layer-integrated pressure differences with face-thickness
         ! divisor.  See OPGF_VARIANT_FV_MOM6 doc above.
      else if (pgf%variant == OPGF_VARIANT_FV_MOM6) then
         if (pgf%reconstruct_for_pressure .and. present(eos) .and. &
             ms%idx_salinity > 0 .and. ms%idx_temperature > 0) then
            ! In-layer PLM/PPM reconstruction: build per-layer Boole
            ! `dpa`/`intz_dpa` from a monotone sub-layer T/S profile, then
            ! reuse the unchanged FV_MOM6 face assembly.
            call compute_fv_mom6_reconstruct_impl(ms%h_layer, &
                                                  ms%tracers(ms%idx_salinity)%hTr, &
                                                  ms%tracers(ms%idx_temperature)%hTr, &
                                                  pgf%b, eos, &
                                                  pgf%recon_S_t%data, pgf%recon_S_b%data, &
                                                  pgf%recon_T_t%data, pgf%recon_T_b%data, &
                                                  pgf%conc_T%data, pgf%conc_S%data, &
                                                  pgf%e_face%data, pgf%pa%data, &
                                                  pgf%intz_dpa%data, &
                                                  pgf%intx_pa%data, pgf%inty_pa%data, &
                                                  pgf%intx_dpa%data, pgf%inty_dpa%data, &
                                                  pgf%dpdx_face%data, pgf%dpdy_face%data, &
                                                  pgf%rho0, pgf%rho_ref, pgf%h_neglect, &
                                                  pgf%gfs_scale, pgf%recon_scheme, &
                                                  ms%p_top, pgf%p_top_in_bc, &
                                                  metrics%idxCu, metrics%idyCv, nx, ny, nz)
         else if (use_insitu_pcm(pgf, ms, eos)) then
            ! Constant-by-layer T/S, density at the in-situ pressure
            ! (MOM6 `int_density_dz_generic_pcm`).  See `insitu_density`.
            call compute_fv_mom6_insitu_pcm_impl(ms%h_layer, &
                                                 ms%tracers(ms%idx_salinity)%hTr, &
                                                 ms%tracers(ms%idx_temperature)%hTr, &
                                                 pgf%b, pgf%conc_T%data, pgf%conc_S%data, &
                                                 pgf%e_face%data, pgf%pa%data, &
                                                 pgf%intz_dpa%data, &
                                                 pgf%intx_pa%data, pgf%inty_pa%data, &
                                                 pgf%intx_dpa%data, pgf%inty_dpa%data, &
                                                 pgf%dpdx_face%data, pgf%dpdy_face%data, &
                                                 pgf%rho0, pgf%rho_ref, pgf%h_neglect, &
                                                 pgf%gfs_scale, pgf%mass_weight, &
                                                 eos%variant, &
                                                 ms%p_top, pgf%p_top_in_bc, &
                                                 metrics%idxCu, metrics%idyCv, nx, ny, nz)
         else
            call compute_fv_mom6_impl(ms%h_layer, ms%rho_layer, pgf%b, &
                                      pgf%e_face%data, pgf%pa%data, &
                                      pgf%intz_dpa%data, &
                                      pgf%intx_pa%data, pgf%inty_pa%data, &
                                      pgf%intx_dpa%data, pgf%inty_dpa%data, &
                                      pgf%dpdx_face%data, pgf%dpdy_face%data, &
                                      pgf%rho0, pgf%rho_ref, pgf%h_neglect, &
                                      pgf%gfs_scale, pgf%mass_weight, &
                                      ms%p_top, pgf%p_top_in_bc, &
                                      metrics%idxCu, metrics%idyCv, nx, ny, nz)
         end if
      else

         ! ---- Pass 1: hydrostatic integration per column ----
         ! Surface boundary condition: p_edge at the top of the
         ! water column is zero (atmospheric absorbed into Boussinesq).
         ! March down: each layer adds rho*g*h_layer to the pressure
         ! at the layer below.  FV_WRIGHT also computes rho_insitu
         ! per layer with a single Picard step.
         if (pgf%variant == OPGF_VARIANT_FV_WRIGHT) then
            if (ms%idx_salinity > 0 .and. ms%idx_temperature > 0) then
               call eos_wright_pgf_column_sweep_impl( &
                  ms%h_layer, &
                  ms%tracers(ms%idx_salinity)%hTr, &
                  ms%tracers(ms%idx_temperature)%hTr, &
                  ms%rho_layer, &
                  ms%p_top, &
                  pgf%p_edge%data, &
                  pgf%rho_insitu%data, &
                  GRAVITY, pgf%rho0, &
                  nx, ny, nz)
            else
               ! No S, T registered: fall through to rho_layer.  Keeps
               ! the test scaffolding (which sets rho_layer directly
               ! without registering tracers) workable.
               do concurrent(j=1:ny, i=1:nx)
                  pgf%p_edge%data(i, j, nz + 1) = 0.0_wp
                  do k = nz, 1, -1
                     pgf%p_edge%data(i, j, k) = pgf%p_edge%data(i, j, k + 1) + &
                                                GRAVITY*ms%rho_layer(i, j, k)*ms%h_layer(i, j, k)
                     pgf%rho_insitu%data(i, j, k) = ms%rho_layer(i, j, k)
                  end do
               end do
            end if
         else
            do concurrent(j=1:ny, i=1:nx)
               pgf%p_edge%data(i, j, nz + 1) = 0.0_wp
               do k = nz, 1, -1
                  pgf%p_edge%data(i, j, k) = pgf%p_edge%data(i, j, k + 1) + &
                                             GRAVITY*ms%rho_layer(i, j, k)*ms%h_layer(i, j, k)
               end do
            end do
         end if

         ! ---- Pass 2: east-face acceleration ----
         ! `idxCu(i,j)` replaces the scalar `inv_dx` (D4) — bit-identical on
         ! uniform Cartesian.
         if (pgf%variant == OPGF_VARIANT_FV_WRIGHT) then
            do concurrent(k=1:nz, j=1:ny, i=2:nx) &
               local(p_centre_left, p_centre_right, rho_face, z_correction)
               p_centre_right = 0.5_wp*(pgf%p_edge%data(i, j, k) + pgf%p_edge%data(i, j, k + 1))
               p_centre_left = 0.5_wp*(pgf%p_edge%data(i - 1, j, k) + pgf%p_edge%data(i - 1, j, k + 1))
               rho_face = 0.5_wp*(pgf%rho_insitu%data(i - 1, j, k) + pgf%rho_insitu%data(i, j, k))
               z_correction = GRAVITY*rho_face* &
                              (pgf%z_centre%data(i, j, k) - pgf%z_centre%data(i - 1, j, k))*metrics%idxCu(i, j)
               pgf%dpdx_face%data(i, j, k) = -inv_rho0*( &
                                             (p_centre_right - p_centre_left)*metrics%idxCu(i, j) + z_correction)
            end do
         else
            ! FV_LITE — the only variant that still reaches this pass
            ! (GPRIME returned, MONT and FV_MOM6 branched above).
            do concurrent(k=1:nz, j=1:ny, i=2:nx) &
               local(p_centre_left, p_centre_right, rho_face, z_correction)
               p_centre_right = 0.5_wp*(pgf%p_edge%data(i, j, k) + pgf%p_edge%data(i, j, k + 1))
               p_centre_left = 0.5_wp*(pgf%p_edge%data(i - 1, j, k) + pgf%p_edge%data(i - 1, j, k + 1))
               rho_face = 0.5_wp*(ms%rho_layer(i - 1, j, k) + ms%rho_layer(i, j, k))
               z_correction = GRAVITY*rho_face* &
                              (pgf%z_centre%data(i, j, k) - pgf%z_centre%data(i - 1, j, k))*metrics%idxCu(i, j)
               pgf%dpdx_face%data(i, j, k) = -inv_rho0*( &
                                             (p_centre_right - p_centre_left)*metrics%idxCu(i, j) + z_correction)
            end do
         end if
         do concurrent(k=1:nz, j=1:ny)
            pgf%dpdx_face%data(1, j, k) = 0.0_wp
            pgf%dpdx_face%data(nx + 1, j, k) = 0.0_wp
         end do

         ! ---- Pass 3: north-face acceleration ----
         ! `idyCv(i,j)` replaces the scalar `inv_dy` (D4).
         if (pgf%variant == OPGF_VARIANT_FV_WRIGHT) then
            do concurrent(k=1:nz, j=2:ny, i=1:nx) &
               local(p_centre_below, p_centre_above, rho_face, z_correction)
               p_centre_above = 0.5_wp*(pgf%p_edge%data(i, j, k) + pgf%p_edge%data(i, j, k + 1))
               p_centre_below = 0.5_wp*(pgf%p_edge%data(i, j - 1, k) + pgf%p_edge%data(i, j - 1, k + 1))
               rho_face = 0.5_wp*(pgf%rho_insitu%data(i, j - 1, k) + pgf%rho_insitu%data(i, j, k))
               z_correction = GRAVITY*rho_face* &
                              (pgf%z_centre%data(i, j, k) - pgf%z_centre%data(i, j - 1, k))*metrics%idyCv(i, j)
               pgf%dpdy_face%data(i, j, k) = -inv_rho0*( &
                                             (p_centre_above - p_centre_below)*metrics%idyCv(i, j) + z_correction)
            end do
         else
            ! FV_LITE — see the Pass-2 comment.
            do concurrent(k=1:nz, j=2:ny, i=1:nx) &
               local(p_centre_below, p_centre_above, rho_face, z_correction)
               p_centre_above = 0.5_wp*(pgf%p_edge%data(i, j, k) + pgf%p_edge%data(i, j, k + 1))
               p_centre_below = 0.5_wp*(pgf%p_edge%data(i, j - 1, k) + pgf%p_edge%data(i, j - 1, k + 1))
               rho_face = 0.5_wp*(ms%rho_layer(i, j - 1, k) + ms%rho_layer(i, j, k))
               z_correction = GRAVITY*rho_face* &
                              (pgf%z_centre%data(i, j, k) - pgf%z_centre%data(i, j - 1, k))*metrics%idyCv(i, j)
               pgf%dpdy_face%data(i, j, k) = -inv_rho0*( &
                                             (p_centre_above - p_centre_below)*metrics%idyCv(i, j) + z_correction)
            end do
         end if
         do concurrent(k=1:nz, i=1:nx)
            pgf%dpdy_face%data(i, 1, k) = 0.0_wp
            pgf%dpdy_face%data(i, ny + 1, k) = 0.0_wp
         end do
      end if

      ! ---- Pass 4: grounded-layer gate (VCOORD_LAGRANGIAN only) ----
      ! Passes 2/3 form a two-point Jacobian: the `Δp_centre` and
      ! `g·ρ_layer·Δz_centre` terms cancel AT REST only while the two abutting
      ! layer centres lie in a common z-interval whose ambient density IS
      ! `ρ_layer`.  Where an isopycnal layer has wedged out against the bed on
      ! one side, the centres are hundreds of metres apart, the interval
      ! between them holds OTHER density classes, and what survives is
      ! `g·(ρ_layer − ρ̄_ambient)·∂z/∂x` — a pressure gradient on a motionless
      ! ocean.  There is no common depth to difference the pressure across
      ! there, so the honest face value is ZERO; the layer is still free to be
      ! re-wetted by continuity's upwind flux and the barotropic correction.
      !
      ! FV_MOM6 is gated by the SAME test.  Its face assembly is not the
      ! two-point Jacobian but the layer-integrated FV-Bouss form, yet the
      ! defect is the geometry, not the quadrature: where the layer occupies
      ! disjoint z-intervals in the two columns the layer-integrated pressure
      ! difference is likewise being taken between depths that share no water
      ! of that density class, and the `1/(h_L + h_R + h_neglect)` divisor
      ! does NOT suppress it — the deep side keeps the thickness up while the
      ! grounded side contributes the whole `e_bot` offset.  Measured on
      ! `seamount_conservative_floor.nml` (form='fv_mom6'): En 2.17e-03 →
      ! 1.39e-26 at day 2.
      !
      ! MONT is gated by the same test for the same reason.  Its face
      ! expression is neither the two-point Jacobian nor the layer-integrated
      ! form, but a grounded layer still puts the two columns' `e_edge` values
      ! hundreds of metres apart, so the `M` recursion stops producing a
      ! horizontally uniform potential at rest and the residual reappears.
      !
      ! GROUNDED, not merely non-overlapping: the face is zeroed only where
      ! the layer is also vanished (`<= nonoverlap_vanish_tol`) on at least
      ! one side.  A layer massive on BOTH sides that merely sits at different
      ! depths (a sigma-seeded stack over a step) carries real mass flux
      ! across the face; its FV PGF is the ordinary steep-coordinate one
      ! (`vcoord_type="sigma"` runs the same geometry), and zeroing it breaks
      ! the PGF-work / PE exchange — see `nonoverlap_vanish_tol`.
      !
      ! A separate guarded pass on purpose: when the gate is off (every vcoord
      ! but LAGRANGIAN) not one extra load is issued ⇒ bit-identical.
      ! Reads `z_centre`, so the driver only ever sets the flag for the four
      ! variants whose allocation gate provides it (MONT / FV_LITE /
      ! FV_WRIGHT, and FV_MOM6 — where `z_centre` is allocated + filled for
      ! this gate alone).
      if (pgf%skip_nonoverlap) then
         vtol = pgf%nonoverlap_vanish_tol
         do concurrent(k=1:nz, j=1:ny, i=2:nx)
            if (min(ms%h_layer(i - 1, j, k), ms%h_layer(i, j, k)) <= vtol) then
               if (min(pgf%z_centre%data(i - 1, j, k) + 0.5_wp*ms%h_layer(i - 1, j, k), &
                       pgf%z_centre%data(i, j, k) + 0.5_wp*ms%h_layer(i, j, k)) <= &
                   max(pgf%z_centre%data(i - 1, j, k) - 0.5_wp*ms%h_layer(i - 1, j, k), &
                       pgf%z_centre%data(i, j, k) - 0.5_wp*ms%h_layer(i, j, k))) then
                  pgf%dpdx_face%data(i, j, k) = 0.0_wp
               end if
            end if
         end do
         do concurrent(k=1:nz, j=2:ny, i=1:nx)
            if (min(ms%h_layer(i, j - 1, k), ms%h_layer(i, j, k)) <= vtol) then
               if (min(pgf%z_centre%data(i, j - 1, k) + 0.5_wp*ms%h_layer(i, j - 1, k), &
                       pgf%z_centre%data(i, j, k) + 0.5_wp*ms%h_layer(i, j, k)) <= &
                   max(pgf%z_centre%data(i, j - 1, k) - 0.5_wp*ms%h_layer(i, j - 1, k), &
                       pgf%z_centre%data(i, j, k) - 0.5_wp*ms%h_layer(i, j, k))) then
                  pgf%dpdy_face%data(i, j, k) = 0.0_wp
               end if
            end if
         end do
      end if
   end subroutine ocean_pressure_force_compute