wright_pcm_dpa_intz Subroutine

   pure subroutine wright_pcm_dpa_intz(t, s, e_top, dz, rho0, rho_ref, dpa, intz_dpa)
      !$acc routine seq
      !! ANALYTIC vertical integral of the Wright (1997) in-situ density
      !! anomaly over one constant-T/S (PCM) layer — MOM6
      !! `int_density_dz_wright` (reduced-range coefficients, MOM6
      !! `EQN_OF_STATE = "WRIGHT"` / `"WRIGHT_RED"`, the set `rdb_eos`
      !! carries).  Replaces the 5-point Boole quadrature of
      !! `boole_dpa_intz_layer` (5 generic-EOS evaluations) with one
      !! polynomial evaluation, one division pair and a short series.
      !!
      !! With `P = p + p0(T,S)`, `L = lambda/alpha0` and the Boussinesq
      !! pressure `p = -g*rho0*z`, Wright's density is
      !!   rho = P/(lambda + alpha0*P) = (1/alpha0)*(1 - L/(P + L)),
      !! so along the layer (T, S fixed)
      !!   int rho dz = dz/alpha0 - (lambda/alpha0^2)/(g*rho0) * ln((1+eps)/(1-eps)),
      !! `eps = (g*rho0*dz/2)/(P_mid + L)` the half-layer pressure change
      !! over the layer-mean `P + L`.  Expanding the log about the
      !! layer midpoint, `ln((1+eps)/(1-eps)) = 2*(eps + eps^3/3 + eps^5/5
      !! + ...)`, the leading term is exactly `dz*rho(P_mid)`, leaving the
      !! remainder `rem = (lambda/alpha0^2)/rho0 * eps^2*(1/3 + eps^2/5 +
      !! eps^4/7 + eps^6/9)`:
      !!   dpa      = g*(rho(P_mid) - rho_ref)*dz - 2*eps*rem
      !!   intz_dpa = 0.5*g*(rho(P_mid) - rho_ref)*dz^2 - dz*(1 + eps)*rem
      !! (`intz_dpa` = the layer integral of the pressure anomaly relative
      !! to its value at the layer TOP — the same moment
      !! `boole_dpa_intz_layer` returns).  The series is truncated after
      !! `eps^8` inside `rem`: `P + L >= ~8e8 Pa` for sea water, so even a
      !! 6000 m layer has `eps < 0.04` and the dropped `eps^10/11` term is
      !! ~1e-15 of `rem` — round-off.  MOM6 states the truncation valid for
      !! `|eps| < 0.34`.
      real(wp), intent(in)  :: t
         !! Layer temperature (degC), constant through the layer.
      real(wp), intent(in)  :: s
         !! Layer salinity (PSU), constant through the layer.
      real(wp), intent(in)  :: e_top
         !! Height of the SHALLOWER interface (m, geopotential, negative
         !! below the datum).
      real(wp), intent(in)  :: dz
         !! Layer thickness (m, >= 0); the layer spans `[e_top - dz, e_top]`.
      real(wp), intent(in)  :: rho0
         !! Boussinesq reference density of the pressure estimate (kg/m^3).
      real(wp), intent(in)  :: rho_ref
         !! Anomaly reference subtracted from the in-situ density (kg/m^3).
      real(wp), intent(out) :: dpa
         !! `g * int (rho - rho_ref) dz` over the layer (Pa).
      real(wp), intent(out) :: intz_dpa
         !! First moment from the top (Pa*m).

      real(wp), parameter :: C1_3 = 1.0_wp/3.0_wp, C1_7 = 1.0_wp/7.0_wp
      real(wp), parameter :: C1_9 = 1.0_wp/9.0_wp
      real(wp) :: al0, p0, lam, gxrho, p_ave, big_p, i_d, half_dp_d, eps, eps2
      real(wp) :: rho_anom, rem

      al0 = WRIGHT_A0 + (WRIGHT_A1*t + WRIGHT_A2*s)
      p0 = WRIGHT_B0 + (WRIGHT_B4*s + t*(WRIGHT_B1 + (t*(WRIGHT_B2 + WRIGHT_B3*t) &
                                                      + WRIGHT_B5*s)))
      lam = WRIGHT_C0 + (WRIGHT_C4*s + t*(WRIGHT_C1 + (t*(WRIGHT_C2 + WRIGHT_C3*t) &
                                                       + WRIGHT_C5*s)))
      ! One division: with D = alpha0*P + lambda (P = p0 + p_ave),
      !   rho(P) = P/D,  1/(P + L) = alpha0/D,
      !   (lambda/alpha0^2)*eps^2 = lambda*(half_dp/D)^2,
      ! algebraically MOM6's `I_al0`/`I_Lzz` form without 1/alpha0.
      gxrho = GRAVITY*rho0
      p_ave = -gxrho*(e_top - 0.5_wp*dz)
      big_p = p0 + p_ave
      i_d = 1.0_wp/(al0*big_p + lam)
      half_dp_d = 0.5_wp*(gxrho*dz)*i_d
      eps = al0*half_dp_d
      eps2 = eps*eps
      rho_anom = big_p*i_d - rho_ref
      rem = (lam/rho0)*(half_dp_d*half_dp_d) &
            *(C1_3 + eps2*(0.2_wp + eps2*(C1_7 + C1_9*eps2)))
      dpa = (GRAVITY*rho_anom)*dz - 2.0_wp*eps*rem
      intz_dpa = 0.5_wp*(GRAVITY*rho_anom)*dz*dz - dz*((1.0_wp + eps)*rem)
   end subroutine wright_pcm_dpa_intz