boole_dpa_intz_layer Subroutine

   pure subroutine boole_dpa_intz_layer(eos, rho0, rho_ref, &
                                        e_top, dz, &
                                        t_t, t_b, t_mean, &
                                        s_t, s_b, s_mean, &
                                        parabolic, dpa, intz_dpa)
      !$acc routine seq
      !! 5-point Boole-quadrature density-anomaly integral over one layer
      !! (Adcroft, Hallberg & Harrison 2008; White, Adcroft & Hallberg
      !! 2009).  Returns the layer-integrated pressure-anomaly increment
      !! `dpa = g * int rho' dz` and its first moment (from the layer
      !! TOP edge inward) `intz_dpa = 0.5 * g * dz^2 * bracket`.
      !!
      !! Sub-layer T/S vary between the top (shallower) edge `_t` and the
      !! bottom (deeper) edge `_b`.  When `parabolic` is .false. (PLM) the
      !! profile is linear in the fractional depth; when .true. (PPM) the
      !! curvature `q6 = 3*(2*q_mean - (q_t + q_b))` is added so the
      !! profile is the in-layer parabola through `_t`, `_b` and the layer
      !! mean.
      !!
      !! The density anomaly at each of the 5 evenly-spaced sub-points
      !! (top -> bottom) is rho'(z) = EOS(T,S,p) - rho_ref with the
      !! Boussinesq hydrostatic pressure estimate p = -g*rho0*z.  rho_ref
      !! is subtracted from the absolute EOS density per sub-point.
      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
         !! Surface-relative height of the SHALLOWER interface (<= 0).
      real(wp), intent(in)  :: dz
         !! Layer thickness (m), dz >= 0.  Marches down from e_top.
      real(wp), intent(in)  :: t_t, t_b, t_mean
         !! Temperature: top edge, bottom edge, layer mean.
      real(wp), intent(in)  :: s_t, s_b, s_mean
         !! Salinity: top edge, bottom edge, layer mean.
      logical, intent(in)  :: parabolic
         !! .true. -> add the PPM curvature (s6/t6) term.
      real(wp), intent(out) :: dpa
         !! g * int rho' dz over the layer (Pa).
      real(wp), intent(out) :: intz_dpa
         !! 0.5 * g * dz^2 * bracket (Pa*m), first moment from the top.

      real(wp) :: gxrho, wt_t, wt_b, t6, s6
      real(wp) :: t5, s5, z5, p5, rho_anom
      real(wp) :: r5(N_BOOLE)
      integer  :: n

      ! The sub-points and weights are written out in place.  The per-EOS
      ! twins of this rule (`boole_dpa_intz_layer_wright`,
      ! `roquet_recon_dpa_intz`) live in `rdb_ocean_pressure_force`, next to
      ! the kernel that calls them, so they inline into it.
      gxrho = GRAVITY*rho0

      ! PPM curvature (zero for PLM).
      t6 = 0.0_wp
      s6 = 0.0_wp
      if (parabolic) then
         t6 = 3.0_wp*(2.0_wp*t_mean - (t_t + t_b))
         s6 = 3.0_wp*(2.0_wp*s_mean - (s_t + s_b))
      end if

      do n = 1, N_BOOLE
         wt_t = 0.25_wp*real(N_BOOLE - n, wp)   ! 1, .75, .5, .25, 0
         wt_b = 1.0_wp - wt_t
         ! Linear blend + parabolic correction.  At wt_t in [0,1] the
         ! parabola through (_t at wt_t=1, _b at wt_t=0, mean) is
         !   q(wt_t) = wt_t*q_t + wt_b*q_b + q6*wt_t*wt_b.
         t5 = wt_t*t_t + wt_b*t_b + t6*wt_t*wt_b
         s5 = wt_t*s_t + wt_b*s_b + s6*wt_t*wt_b
         z5 = e_top - 0.25_wp*real(n - 1, wp)*dz   ! marches DOWN from top
         p5 = -gxrho*z5
         r5(n) = eos_density_point(eos, t5, s5, p5) - rho_ref
      end do

      rho_anom = (1.0_wp/90.0_wp)*(7.0_wp*(r5(1) + r5(5)) &
                                   + 32.0_wp*(r5(2) + r5(4)) + 12.0_wp*r5(3))
      dpa = GRAVITY*dz*rho_anom
      intz_dpa = 0.5_wp*GRAVITY*dz*dz*(rho_anom &
                                       - (1.0_wp/90.0_wp)*(16.0_wp*(r5(4) - r5(2)) + 7.0_wp*(r5(5) - r5(1))))
   end subroutine boole_dpa_intz_layer