wright_pcm_dpa_face Subroutine

   pure subroutine wright_pcm_dpa_face(e_top_l, e_top_r, dz_l, dz_r, &
                                       t_l, t_r, s_l, s_r, dpa_l, dpa_r, &
                                       hwt_ll, hwt_lr, hwt_rr, hwt_rl, &
                                       rho0, rho_ref, dpa_face)
      !$acc routine seq
      !! Cross-face 5-point Boole quadrature of the layer `dpa` for a PCM
      !! column pair with the ANALYTIC Wright vertical integral at every
      !! lateral sub-column (MOM6 `int_density_dz_wright`, `intx_dpa` /
      !! `inty_dpa`).  The Wright twin of `boole_dpa_face_pcm`, with the
      !! same sub-column definition: the end points are the columns' own
      !! `dpa_l` / `dpa_r`; the three interior sub-columns interpolate the
      !! interface height and thickness linearly in the cross-face
      !! fraction and T/S with the near-bottom mass-weighted fractions
      !! `wtT_L = wl*hwt_ll + wr*hwt_rl`, `wtT_R = wl*hwt_lr + wr*hwt_rr`.
      !!
      !! (MOM6 interpolates the Wright polynomial TERMS `alpha0, p0,
      !! lambda` rather than T and S.  Interpolating T/S keeps the
      !! sub-column identical to the Boole path's, so switching the vertical
      !! rule is the only change; the two differ by the polynomials'
      !! curvature, `O(wl*wr*dT^2)`, which is zero on a resting column.)
      real(wp), intent(in)  :: e_top_l, e_top_r
         !! Shallower-interface heights in the left / right column (m).
      real(wp), intent(in)  :: dz_l, dz_r
         !! Layer thicknesses in the left / right column (m, >= 0).
      real(wp), intent(in)  :: t_l, t_r, s_l, s_r
         !! Layer-mean temperature / salinity in the left / right column.
      real(wp), intent(in)  :: dpa_l, dpa_r
         !! The columns' own `g * int rho' dz` over the layer (Pa).
      real(wp), intent(in)  :: hwt_ll, hwt_lr, hwt_rr, hwt_rl
         !! MOM6 `hWt_LL/LR/RR/RL` mass-weighting fractions.
      real(wp), intent(in)  :: rho0, rho_ref
         !! Boussinesq pressure-estimate density / anomaly reference.
      real(wp), intent(out) :: dpa_face
         !! Along-face mean of `g * int rho' dz` over the layer (Pa).

      real(wp) :: wr, wl, wtt_l, wtt_r, dpa_m, intz_m, acc
      integer  :: m

      acc = 7.0_wp*(dpa_l + dpa_r)
      ! (`GCC$ unroll`: gfortran unrolls the three sub-columns completely
      ! before vectorising, so the face loops over this vectorise on the CPU;
      ! nvfortran keeps the loop, which costs fewer registers on the device.)
      !GCC$ unroll 3
      do m = 2, 4
         wr = 0.25_wp*real(m - 1, wp)   ! 0 at the left column .. 1 at the right
         wl = 1.0_wp - wr
         wtt_l = wl*hwt_ll + wr*hwt_rl
         wtt_r = wl*hwt_lr + wr*hwt_rr
         call wright_pcm_dpa_intz(wtt_l*t_l + wtt_r*t_r, wtt_l*s_l + wtt_r*s_r, &
                                  wl*e_top_l + wr*e_top_r, wl*dz_l + wr*dz_r, &
                                  rho0, rho_ref, dpa_m, intz_m)
         if (m == 3) then
            acc = acc + 12.0_wp*dpa_m
         else
            acc = acc + 32.0_wp*dpa_m
         end if
      end do
      dpa_face = acc/90.0_wp
   end subroutine wright_pcm_dpa_face