compute_fv_mom6_impl Subroutine

   pure subroutine compute_fv_mom6_impl(h_layer, rho_layer, b, &
                                        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, &
                                        p_top, p_top_in_bc, &
                                        idxCu, idyCv, nx, ny, nz)
      !! Faithful port of MOM6's `PressureForce_FV_Bouss` per-layer PGF
      !! for the Boussinesq + per-layer Rlay path.
      !!
      !! Pass layout (matches MOM6's `PressureForce_FV_Bouss`):
      !!
      !!   Pass 1 (per column):
      !!     e_face(i, j, k_face) — interface heights, positive up.
      !!       e_face(1) = -b (bed); e_face(k+1) = e_face(k) + h_layer(k);
      !!       e_face(nz+1) = -b + sum(h_layer) = η (free surface).
      !!     pa(i, j, nz+1) = rho_ref · g · η  (surface BC), plus
      !!       `p_top(i, j)` when `p_top_in_bc` (see the theorem below).
      !!     pa(i, j, k) = pa(i, j, k+1) + (rho_layer(k) − rho_ref) · g · h(k)
      !!       (marching down).
      !!     intz_dpa(i, j, k) = 0.5 · (rho_layer(k) − rho_ref) · g · h(k)²
      !!       (mid-point rule).
      !!
      !!   Pass 2 (per face, march down):
      !!     intx_pa(i, j, nz+1) = 0.5 · (pa(i-1, ., nz+1) + pa(i, ., nz+1))
      !!       (surface BC).
      !!     intx_dpa(i, j, k) = 0.5 · g · ((rho(i-1, k) − rho_ref) · h(i-1, k)
      !!                                    + (rho(i, k) − rho_ref) · h(i, k))
      !!       (generalised to per-cell rho).
      !!     intx_pa(i, j, k) = intx_pa(i, j, k+1) + intx_dpa(i, j, k)
      !!       (face pressure recurrence).
      !!     Symmetric on v-face.
      !!
      !!     When `mass_weight = .true.`, at hydrostatically-inconsistent
      !!     unequal-depth faces (`hWght > 0`,
      !!     `hWght = max(0, e_bed_R − e_top_L,k, e_bed_L − e_top_R,k)`)
      !!     the layer density entering `dpa_L`/`dpa_R` is replaced by
      !!     the MOM6 `hWt_LL/LR/RR/RL` blend biased toward the
      !!     thinner column, AND the integral uses the face-interpolated
      !!     thickness `dz = 0.5·(h_L + h_R)` for BOTH samples (MOM6
      !!     `dz_x · rho_anom`).  The interpolated
      !!     thickness is load-bearing: the plain per-cell form
      !!     `0.5·g·(ρ_L'·h_L + ρ_R'·h_R)` is invariant under the
      !!     (Σρh-conserving) hWt blend, so the shelf-break cancellation
      !!     only appears when one `dz` multiplies both blended
      !!     densities.  `hWght = 0` (aligned / equal-depth) ⇒ the exact
      !!     per-cell layer-midpoint average ⇒ bit-identical.
      !!
      !!   Pass 3 (PFu/PFv assembly):
      !!     numer = ((pa(L, k+1) · h(L, k) + intz_dpa(L, k))
      !!              − (pa(R, k+1) · h(R, k) + intz_dpa(R, k)))
      !!           + (h(R, k) − h(L, k)) · intx_pa(face, k+1)
      !!           − (e_face(R, k) − e_face(L, k)) · intx_dpa(face, k)
      !!     denom = h(L, k) + h(R, k) + h_neglect
      !!     PFu(face, k) = numer · (2 · I_Rho0 · IdxCu) / denom
      !!
      !! Convention map (MOM6 → Roundabout):
      !!   MOM6 K (top of layer k_mom6)        → ours k+1 (top of layer k)
      !!   MOM6 K+1 (bottom of layer k_mom6)   → ours k   (bottom of layer k)
      !!   MOM6 i (left of u-face I)           → ours i-1 (west of u-face i)
      !!   MOM6 i+1 (right of u-face I)        → ours i   (east of u-face i)
      !!   MOM6 k_mom6 = 1 (surface layer)     → ours k = nz
      !!   MOM6 k_mom6 = nz_mom6 (bed layer)   → ours k = 1
      !!
      !! Wall faces (face index 1 and N+1) get zero by convention — the
      !! BT-substep / slow continuity already enforces u=0 there.
      !!
      !! ## Theorem — a depth-uniform top load is baroclinically inert here
      !!
      !! Perturb the TOP boundary condition only: `pa(·,nz+1) → pa(·,nz+1)
      !! + δp` with `δp(i,j)` independent of `k`. Every `dpa`, `intz_dpa`
      !! and `intx_dpa` is unchanged, and the Pass-2 recurrence shifts
      !! `intx_pa(K) → intx_pa(K) + ½(δp_L + δp_R)` for EVERY `K`. The
      !! Pass-3 numerator therefore moves by
      !!
      !!   δnumer = δp_L·h_L − δp_R·h_R + (h_R − h_L)·½(δp_L + δp_R)
      !!          = ½(h_L + h_R)·(δp_L − δp_R)
      !!   δPFu(k) = −(1/ρ₀)·(δp_R − δp_L)·IdxCu · (h_L+h_R)/(h_L+h_R+h_neglect)
      !!
      !! — i.e. exactly `−(1/ρ₀)·∂δp/∂x`, **the same in every layer**, up
      !! to the `h_neglect` divisor (a relative `h_n/h_k ≈ 1e-10` for a
      !! metre-thick layer, `≈7e-7` for one at `H_VANISHED`).
      !!
      !! Consequence for the SPLIT solver (`&ocean_bt_nml bc_pgf_forcing`,
      !! default, MOM6 `BT_force`): the depth mean of the layer PGF FORCES
      !! the barotropic substep, so a depth-uniform `δPFu` reaches the
      !! barotropic mode.  The part of `p_top` that is the atmospheric /
      !! anomaly load `sf%p_surf` is ALSO on the `eta_forcing` seam, and
      !! `set_fast_forcing_eta_pf` sheds `g·∇η_ib` from the forcing so it
      !! is counted once; the static ice load `p_ice_ref` cancels inside
      !! `pa(nz+1)` against the datum-shifted `η_geo`, and what survives
      !! of it is the physical reference-density shortfall.  (The legacy
      !! split, `bc_pgf_forcing = .false.`, subtracted the whole depth
      !! mean, so there a depth-uniform `δPFu` cancelled identically.)
      !! The UNSPLIT driver has no seam, so there this term is the load's
      !! only path into the momentum.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in)    :: h_layer(nx, ny, nz)
      real(wp), intent(in)    :: rho_layer(nx, ny, nz)
      real(wp), intent(in)    :: b(nx, ny)
      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
      real(wp), intent(in)    :: p_top(nx, ny)
         !! Top-of-column pressure (Pa, `>= 0`), `multilayer_state_t%p_top`.
         !! Consulted only when `p_top_in_bc`; the zero array otherwise.
      logical, intent(in)    :: p_top_in_bc
         !! Add `p_top` to the Pass-1 surface BC. `.false.` ⇒ the
         !! assignment is character-for-character the pre-knob one ⇒
         !! bit-identical (same branch-on-a-scalar-knob shape as
         !! `mass_weight` in Pass 2).
      real(wp), intent(in)    :: idxCu(nx + 1, ny), idyCv(nx, ny + 1)

      integer  :: i, j, k
      real(wp) :: inv_rho0, eta, rho_anom, dpa_kk
      real(wp) :: dpa_L, dpa_R, 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
      real(wp) :: hwght, hwl, hwr, idenom_hw, hwt_ll, hwt_lr, hwt_rr, hwt_rl
      real(wp) :: rho_face_l, rho_face_r

      inv_rho0 = 1.0_wp/rho0

      ! ---- Pass 1: per-column build of e_face, pa, intz_dpa ----
      do concurrent(j=1:ny, i=1:nx) local(k, eta, rho_anom, dpa_kk)
         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
         eta = e_face(i, j, nz + 1)
         if (p_top_in_bc) then
            pa(i, j, nz + 1) = rho_ref*GRAVITY*eta + p_top(i, j)
         else
            pa(i, j, nz + 1) = rho_ref*GRAVITY*eta
         end if
         do k = nz, 1, -1
            rho_anom = rho_layer(i, j, k) - rho_ref
            dpa_kk = rho_anom*GRAVITY*h_layer(i, j, k)
            pa(i, j, k) = pa(i, j, k + 1) + dpa_kk
            intz_dpa(i, j, k) = 0.5_wp*dpa_kk*h_layer(i, j, k)
         end do
      end do

      ! ---- Pass 2a: u-face horizontal integrals ----
      do concurrent(j=1:ny, i=2:nx) local(k, dpa_L, dpa_R, hwght, hwl, hwr, &
                                          idenom_hw, hwt_ll, hwt_lr, hwt_rr, hwt_rl, &
                                          rho_face_l, rho_face_r)
         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
            ! hWght: hydrostatic-inconsistency measure at this u-face for
            ! layer k.  Bed height = e_face(.,1); layer-k top = e_face(.,k+1).
            ! MOM6 form: max(0, bed_R - top_L, bed_L - top_R).
            hwght = 0.0_wp
            if (mass_weight) then
               hwght = max(0.0_wp, &
                           e_face(i, j, 1) - e_face(i - 1, j, k + 1), &
                           e_face(i - 1, j, 1) - e_face(i, j, k + 1))
            end if
            if (hwght > 0.0_wp) then
               ! Hydrostatically-inconsistent face: blend the layer
               ! density toward the thinner column (MOM6 hWt_*) and
               ! integrate with the face-interpolated thickness
               ! dz = 0.5·(h_L + h_R).  The interpolated thickness is
               ! what makes the blend non-trivial — the plain per-cell
               ! form `0.5·g·(ρ_L'·h_L + ρ_R'·h_R)` is invariant under
               ! the (Σρh-conserving) hWt blend, so the cancellation
               ! only appears when the SAME dz multiplies both samples
               ! (MOM6 `dz_x · rho_anom`).
               hwl = h_layer(i - 1, j, k) + h_neglect
               hwr = h_layer(i, j, k) + h_neglect
               hwght = hwght*((hwl - hwr)/(hwl + hwr))**2
               idenom_hw = 1.0_wp/(hwght*(hwr + hwl) + hwl*hwr)
               hwt_ll = (hwght*hwl + hwr*hwl)*idenom_hw
               hwt_lr = (hwght*hwr)*idenom_hw
               hwt_rr = (hwght*hwr + hwr*hwl)*idenom_hw
               hwt_rl = (hwght*hwl)*idenom_hw
               rho_face_l = hwt_ll*rho_layer(i - 1, j, k) + hwt_lr*rho_layer(i, j, k)
               rho_face_r = hwt_rl*rho_layer(i - 1, j, k) + hwt_rr*rho_layer(i, j, k)
               dpa_L = (rho_face_l - rho_ref)*GRAVITY* &
                       (0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k)))
               dpa_R = (rho_face_r - rho_ref)*GRAVITY* &
                       (0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k)))
            else
               ! Aligned / equal-depth face: exact layer-midpoint
               ! average (bit-identical to the pre-knob FV_MOM6 form).
               dpa_L = (rho_layer(i - 1, j, k) - rho_ref)*GRAVITY*h_layer(i - 1, j, k)
               dpa_R = (rho_layer(i, j, k) - rho_ref)*GRAVITY*h_layer(i, j, k)
            end if
            intx_dpa(i, j, k) = 0.5_wp*(dpa_L + dpa_R)
            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 ----
      do concurrent(j=2:ny, i=1:nx) local(k, dpa_L, dpa_R, hwght, hwl, hwr, &
                                          idenom_hw, hwt_ll, hwt_lr, hwt_rr, hwt_rl, &
                                          rho_face_l, rho_face_r)
         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
            hwght = 0.0_wp
            if (mass_weight) then
               hwght = max(0.0_wp, &
                           e_face(i, j, 1) - e_face(i, j - 1, k + 1), &
                           e_face(i, j - 1, 1) - e_face(i, j, k + 1))
            end if
            if (hwght > 0.0_wp) then
               ! Hydrostatically-inconsistent face: hWt density blend +
               ! face-interpolated thickness (see Pass 2a comment).
               hwl = h_layer(i, j - 1, k) + h_neglect
               hwr = h_layer(i, j, k) + h_neglect
               hwght = hwght*((hwl - hwr)/(hwl + hwr))**2
               idenom_hw = 1.0_wp/(hwght*(hwr + hwl) + hwl*hwr)
               hwt_ll = (hwght*hwl + hwr*hwl)*idenom_hw
               hwt_lr = (hwght*hwr)*idenom_hw
               hwt_rr = (hwght*hwr + hwr*hwl)*idenom_hw
               hwt_rl = (hwght*hwl)*idenom_hw
               rho_face_l = hwt_ll*rho_layer(i, j - 1, k) + hwt_lr*rho_layer(i, j, k)
               rho_face_r = hwt_rl*rho_layer(i, j - 1, k) + hwt_rr*rho_layer(i, j, k)
               dpa_L = (rho_face_l - rho_ref)*GRAVITY* &
                       (0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k)))
               dpa_R = (rho_face_r - rho_ref)*GRAVITY* &
                       (0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k)))
            else
               ! Aligned / equal-depth face: exact layer-midpoint average.
               dpa_L = (rho_layer(i, j - 1, k) - rho_ref)*GRAVITY*h_layer(i, j - 1, k)
               dpa_R = (rho_layer(i, j, k) - rho_ref)*GRAVITY*h_layer(i, j, k)
            end if
            inty_dpa(i, j, k) = 0.5_wp*(dpa_L + dpa_R)
            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 ----
      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) ----
      ! Subtracts `(1 - gfs_scale)·(g/ρ₀)·ρ_surf·∇η` from every layer's
      ! PGF so the slow tendency only carries `gfs_scale · g · ∇η`.
      ! The BT substep evolves η with `g_bt = gfs_scale · GRAVITY`
      ! (set by the driver); combined they reproduce the physical
      ! surface-gravity coupling at the chosen reduced value.
      ! Depth-independent → BT mass-flux invariant intact.  Mirrors
      ! MOM6's Boussinesq non-EOS branch
      ! (`rho_surf = Rlay(top)`).  No-op when `gfs_scale = 1`
      ! (modulo the small tolerance below) so the default
      ! configuration is bit-identical to the pre-knob FV_MOM6 port.
      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*(rho_layer(i, j, nz)*e_face(i, j, nz + 1) &
                               - rho_layer(i - 1, j, nz)*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*(rho_layer(i, j, nz)*e_face(i, j, nz + 1) &
                               - rho_layer(i, j - 1, nz)*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_impl