boole_dpa_face Subroutine

   pure subroutine boole_dpa_face(eos, rho0, rho_ref, &
                                  e_top_l, e_top_r, dz_l, dz_r, &
                                  t_t_l, t_b_l, t_m_l, t_t_r, t_b_r, t_m_r, &
                                  s_t_l, s_b_l, s_m_l, s_t_r, s_b_r, s_m_r, &
                                  dpa_l, dpa_r, parabolic, dpa_face)
      !$acc routine seq
      !! HORIZONTAL (cross-face) Boole quadrature of the layer pressure
      !! increment `dpa = g * int rho' dz` — the face integral the FV
      !! pressure-gradient contour needs (Adcroft, Hallberg & Harrison
      !! 2008 §3; Yung, Hallberg, Adcroft & Morrison 2026 §2.4).
      !!
      !! WHY A QUADRATURE AND NOT A MEAN.  The FV assembly evaluates the
      !! top/bottom edges of the control volume as `Delta_e * pbar`, where
      !! `pbar` is the mean pressure ALONG that edge; `pbar` is marched
      !! down from the surface by adding this routine's result layer by
      !! layer.  `Delta_e * pbar` is the exact `int p dz` along the edge
      !! only when `pbar` is the true along-face mean.  Replacing it with
      !! the two-column average `0.5*(dpa_L + dpa_R)` is a TRAPEZOID: it is
      !! exact only if `p` is linear in x along the edge.  With a tilted
      !! interface, `z` is linear in x but `p` is QUADRATIC in z under a
      !! linear stratification, so the trapezoid leaves a curvature
      !! residual `g*(-drho/dz)*Delta_e^2/12` at every interface — the
      !! terrain-following "pressure gradient error of the second kind"
      !! (Haney 1991; Mellor, Ezer & Oey 1994), reported for the sloping
      !! ice-shelf surface as the LINEAR pressure reconstruction by Yung
      !! et al. (2026) §3.1 and cured there by the same device.
      !!
      !! THE FIX.  Sample five evenly-spaced sub-columns across the face.
      !! At fraction `w` from the left column, interpolate the interface
      !! height, the thickness, and the T/S edge + mean triples linearly,
      !! then run the same in-layer vertical Boole quadrature.  Both
      !! interpolations are in the SAME parameter, so a sub-column's
      !! profile is the true profile at the interpolated depth: for T/S
      !! linear in z the sub-column reconstruction is exact and `dpa(w)`
      !! is a quadratic in `w`, which the 5-point closed Newton-Cotes rule
      !! integrates exactly.  The whole PGF then vanishes to round-off on
      !! a resting linear-EOS/linear-stratification column under ANY
      !! layer geometry — the algorithm's defining property.
      !!
      !! THE END POINTS ARE THE COLUMNS' OWN INTEGRALS.  At `w = 0` and
      !! `w = 1` the sub-column IS the left / right column, so its `dpa` is
      !! the one Pass 1 already integrated; the caller passes it in
      !! (`dpa_l` / `dpa_r`, recovered from the `pa` stack as MOM6
      !! `int_density_dz_generic_plm` does) and only the three interior
      !! sub-columns are integrated here.
      !!
      !! COST: 3 sub-columns x 5 sub-points = 15 EOS evaluations per face
      !! per layer (25 before the end points were reused), against 0 for the
      !! trapezoid.  `reconstruct_for_pressure` is opt-in and already the
      !! expensive branch.
      type(eos_t), intent(in) :: eos
      real(wp), intent(in)  :: rho0
         !! Boussinesq reference density used in the pressure estimate.
      real(wp), intent(in)  :: rho_ref
         !! Anomaly reference subtracted from the EOS density.
      real(wp), intent(in)  :: e_top_l, e_top_r
         !! Shallower-interface heights of the layer in the LEFT and
         !! RIGHT columns (m, surface-relative, <= 0).
      real(wp), intent(in)  :: dz_l, dz_r
         !! Layer thicknesses in the left / right columns (m, >= 0).
      real(wp), intent(in)  :: t_t_l, t_b_l, t_m_l
         !! Left column temperature: top edge, bottom edge, layer mean.
      real(wp), intent(in)  :: t_t_r, t_b_r, t_m_r
         !! Right column temperature triple.
      real(wp), intent(in)  :: s_t_l, s_b_l, s_m_l
         !! Left column salinity triple.
      real(wp), intent(in)  :: s_t_r, s_b_r, s_m_r
         !! Right column salinity triple.
      real(wp), intent(in)  :: dpa_l, dpa_r
         !! The left / right columns' own `g * int rho' dz` over the layer
         !! (Pa) -- the `w = 0` / `w = 1` end points of the rule.
      logical, intent(in)  :: parabolic
         !! .true. -> the sub-column profiles carry the PPM curvature.
      real(wp), intent(out) :: dpa_face
         !! Along-face mean of `g * int rho' dz` over the layer (Pa).

      real(wp) :: wr, wl, dpa_m, intz_m, acc
      integer  :: m
      real(wp), parameter :: BOOLE_W(N_BOOLE) = &
                             [7.0_wp, 32.0_wp, 12.0_wp, 32.0_wp, 7.0_wp]

      acc = BOOLE_W(1)*dpa_l + BOOLE_W(N_BOOLE)*dpa_r
      do m = 2, N_BOOLE - 1
         wr = 0.25_wp*real(m - 1, wp)   ! 0 at the left column .. 1 at the right
         wl = 1.0_wp - wr
         call boole_dpa_intz_layer(eos, rho0, rho_ref, &
                                   wl*e_top_l + wr*e_top_r, &
                                   wl*dz_l + wr*dz_r, &
                                   wl*t_t_l + wr*t_t_r, &
                                   wl*t_b_l + wr*t_b_r, &
                                   wl*t_m_l + wr*t_m_r, &
                                   wl*s_t_l + wr*s_t_r, &
                                   wl*s_b_l + wr*s_b_r, &
                                   wl*s_m_l + wr*s_m_r, &
                                   parabolic, dpa_m, intz_m)
         acc = acc + BOOLE_W(m)*dpa_m
      end do
      dpa_face = acc/90.0_wp
   end subroutine boole_dpa_face