MOM_density_integrals.F90

!> Provides integrals of density
module mom_density_integrals

! This file is part of MOM6. See LICENSE.md for the license.

use mom_eos,              only : eos_type
use mom_eos,              only : eos_quadrature
use mom_eos,              only : analytic_int_density_dz
use mom_eos,              only : analytic_int_specific_vol_dp
use mom_eos,              only : calculate_density
use mom_eos,              only : calculate_spec_vol
use mom_eos,              only : calculate_specific_vol_derivs
use mom_error_handler,    only : mom_error, fatal, warning, mom_mesg
use mom_file_parser,      only : get_param, log_version, param_file_type
use mom_hor_index,        only : hor_index_type
use mom_string_functions, only : uppercase
use mom_variables,        only : thermo_var_ptrs
use mom_unit_scaling,     only : unit_scale_type
use mom_verticalgrid,     only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public int_density_dz
public int_density_dz_generic_pcm
public int_density_dz_generic_plm
public int_density_dz_generic_ppm
public int_specific_vol_dp
public int_spec_vol_dp_generic_pcm
public int_spec_vol_dp_generic_plm
public find_depth_of_pressure_in_cell

contains

!> Calls the appropriate subroutine to calculate analytical and nearly-analytical
!! integrals in z across layers of pressure anomalies, which are
!! required for calculating the finite-volume form pressure accelerations in a
!! Boussinesq model.
subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, &
                          intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hi  !< Ocean horizontal index structures for the arrays
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: t   !< Potential temperature referenced to the surface [degC]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: s   !< Salinity [ppt]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: z_t !< Height at the top of the layer in depth units [Z ~> m]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m]
  real,                 intent(in)  :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is
                                           !! subtracted out to reduce the magnitude of each of the
                                           !! integrals.
  real,                 intent(in)  :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used
                                           !! to calculate the pressure (as p~=-z*rho_0*G_e)
                                           !! used in the equation of state.
  real,                 intent(in)  :: g_e !< The Earth's gravitational acceleration
                                           !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2]
  type(eos_type),       pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in) :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(HI),SZJ_(HI)), &
                      intent(inout) :: dpa !< The change in the pressure anomaly
                                           !! across the layer [R L2 T-2 ~> Pa] or [Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
            optional, intent(inout) :: intz_dpa !< The integral through the thickness of the
                                           !! layer of the pressure anomaly relative to the
                                           !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]
  real, dimension(SZIB_(HI),SZJ_(HI)), &
            optional, intent(inout) :: intx_dpa !< The integral in x of the difference between
                                          !! the pressure anomaly at the top and bottom of the
                                          !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJB_(HI)), &
            optional, intent(inout) :: inty_dpa !< The integral in y of the difference between
                                          !! the pressure anomaly at the top and bottom of the
                                          !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(in)  :: bathyt !< The depth of the bathymetry [Z ~> m]
  real,       optional, intent(in)  :: dz_neglect !< A minuscule thickness change [Z ~> m]
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting to
                                           !! interpolate T/S for top and bottom integrals.

  if (eos_quadrature(eos)) then
    call int_density_dz_generic_pcm(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, eos, us, dpa, &
                                    intz_dpa, intx_dpa, inty_dpa, bathyt, dz_neglect, usemasswghtinterp)
  else
    call analytic_int_density_dz(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, eos, dpa, &
                                 intz_dpa, intx_dpa, inty_dpa, bathyt, dz_neglect, usemasswghtinterp)
  endif

end subroutine int_density_dz


!> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which
!! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model.
subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
                                      EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
                                      bathyT, dz_neglect, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hi  !< Horizontal index type for input variables.
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: t  !< Potential temperature of the layer [degC]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: s  !< Salinity of the layer [ppt]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: z_t !< Height at the top of the layer in depth units [Z ~> m]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m]
  real,                 intent(in)  :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is
                                          !! subtracted out to reduce the magnitude
                                          !! of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used
                                          !! to calculate the pressure (as p~=-z*rho_0*G_e)
                                          !! used in the equation of state.
  real,                 intent(in)  :: g_e !< The Earth's gravitational acceleration
                                          !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2]
  type(eos_type),       pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
  real, dimension(SZI_(HI),SZJ_(HI)), &
                      intent(inout) :: dpa !< The change in the pressure anomaly
                                          !! across the layer [R L2 T-2 ~> Pa] or [Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
            optional, intent(inout) :: intz_dpa !< The integral through the thickness of the
                                          !! layer of the pressure anomaly relative to the
                                          !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]
  real, dimension(SZIB_(HI),SZJ_(HI)), &
            optional, intent(inout) :: intx_dpa !< The integral in x of the difference between
                                          !! the pressure anomaly at the top and bottom of the
                                          !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJB_(HI)), &
            optional, intent(inout) :: inty_dpa !< The integral in y of the difference between
                                          !! the pressure anomaly at the top and bottom of the
                                          !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(in)  :: bathyt !< The depth of the bathymetry [Z ~> m]
  real,       optional, intent(in)  :: dz_neglect !< A minuscule thickness change [Z ~> m]
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting to
                                          !! interpolate T/S for top and bottom integrals.
  ! Local variables
  real :: t5(5), s5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt]
  real :: p5(5)      ! Pressures at five quadrature points, never rescaled from Pa [Pa]
  real :: r5(5)      ! Densities at five quadrature points [R ~> kg m-3] or [kg m-3]
  real :: rho_anom   ! The depth averaged density anomaly [R ~> kg m-3] or [kg m-3]
  real :: w_left, w_right ! Left and right weights [nondim]
  real, parameter :: c1_90 = 1.0/90.0  ! Rational constants.
  real :: gxrho      ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
  real :: i_rho      ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]
  real :: rho_scale  ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]
  real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]
  real :: dz         ! The layer thickness [Z ~> m]
  real :: hwght      ! A pressure-thickness below topography [Z ~> m]
  real :: hl, hr     ! Pressure-thicknesses of the columns to the left and right [Z ~> m]
  real :: idenom     ! The inverse of the denominator in the weights [Z-2 ~> m-2]
  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]
  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]
  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]
  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]
  real :: intz(5)    ! The gravitational acceleration times the integrals of density
                     ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]
  logical :: do_massweight ! Indicates whether to do mass weighting.
  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n

  ! These array bounds work for the indexing convention of the input arrays, but
  ! on the computational domain defined for the output arrays.
  isq = hi%IscB ; ieq = hi%IecB
  jsq = hi%JscB ; jeq = hi%JecB
  is = hi%isc ; ie = hi%iec
  js = hi%jsc ; je = hi%jec

  rho_scale = us%kg_m3_to_R
  gxrho = us%RL2_T2_to_Pa * g_e * rho_0
  rho_ref_mks = rho_ref * us%R_to_kg_m3
  i_rho = 1.0 / rho_0

  do_massweight = .false.
  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
    do_massweight = .true.
    if (.not.present(bathyt)) call mom_error(fatal, "int_density_dz_generic: "//&
        "bathyT must be present if useMassWghtInterp is present and true.")
    if (.not.present(dz_neglect)) call mom_error(fatal, "int_density_dz_generic: "//&
        "dz_neglect must be present if useMassWghtInterp is present and true.")
  endif ; endif

  do j=jsq,jeq+1 ; do i=isq,ieq+1
    dz = z_t(i,j) - z_b(i,j)
    do n=1,5
      t5(n) = t(i,j) ; s5(n) = s(i,j)
      p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
    enddo
    if (rho_scale /= 1.0) then
      call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
    else
      call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks)
    endif

    ! Use Boole's rule to estimate the pressure anomaly change.
    rho_anom = c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))
    dpa(i,j) = g_e*dz*rho_anom
    ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
    ! the pressure anomaly.
    if (present(intz_dpa)) intz_dpa(i,j) = 0.5*g_e*dz**2 * &
          (rho_anom - c1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
  enddo ; enddo

  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
    if (hwght > 0.) then
      hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
      hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif

    intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
    do m=2,4
      ! T, S, and z are interpolated in the horizontal.  The z interpolation
      ! is linear, but for T and S it may be thickness weighted.
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
      dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
      p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i+1,j))
      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo
      if (rho_scale /= 1.0) then
        call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks)
      endif

    ! Use Boole's rule to estimate the pressure anomaly change.
      intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
    enddo
    ! Use Boole's rule to integrate the bottom pressure anomaly values in x.
    intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                           12.0*intz(3))
  enddo ; enddo ; endif

  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
    if (hwght > 0.) then
      hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
      hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif

    intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
    do m=2,4
      ! T, S, and z are interpolated in the horizontal.  The z interpolation
      ! is linear, but for T and S it may be thickness weighted.
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
      dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
      p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i,j+1))
      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1)
        p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo
      if (rho_scale /= 1.0) then
        call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks)
      endif

    ! Use Boole's rule to estimate the pressure anomaly change.
      intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
    enddo
    ! Use Boole's rule to integrate the values.
    inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                                     12.0*intz(3))
  enddo ; enddo ; endif
end subroutine int_density_dz_generic_pcm


!> Compute pressure gradient force integrals by quadrature for the case where
!! T and S are linear profiles.
subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
                                      rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, dpa, &
                                      intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp)
  integer,              intent(in)  :: k   !< Layer index to calculate integrals for
  type(hor_index_type), intent(in)  :: hi  !< Ocean horizontal index structures for the input arrays
  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
  type(thermo_var_ptrs), intent(in) :: tv  !< Thermodynamic variables
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: t_t !< Potential temperature at the cell top [degC]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: t_b !< Potential temperature at the cell bottom [degC]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: s_t !< Salinity at the cell top [ppt]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: s_b !< Salinity at the cell bottom [ppt]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &
                        intent(in)  :: e   !< Height of interfaces [Z ~> m]
  real,                 intent(in)  :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is subtracted
                                           !! out to reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used to calculate
                                           !! the pressure (as p~=-z*rho_0*G_e) used in the equation of state.
  real,                 intent(in)  :: g_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
  real,                 intent(in)  :: dz_subroundoff !< A minuscule thickness change [Z ~> m]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: bathyt !< The depth of the bathymetry [Z ~> m]
  type(eos_type),       pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [R L2 Z T-2 ~> Pa Z]
  real, dimension(SZIB_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the x grid spacing [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJB_(HI)), &
              optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the y grid spacing [R L2 T-2 ~> Pa]
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting to
                                           !! interpolate T/S for top and bottom integrals.

! This subroutine calculates (by numerical quadrature) integrals of
! pressure anomalies across layers, which are required for calculating the
! finite-volume form pressure accelerations in a Boussinesq model.  The one
! potentially dodgy assumption here is that rho_0 is used both in the denominator
! of the accelerations, and in the pressure used to calculated density (the
! latter being -z*rho_0*G_e).  These two uses could be separated if need be.
!
! It is assumed that the salinity and temperature profiles are linear in the
! vertical. The top and bottom values within each layer are provided and
! a linear interpolation is used to compute intermediate values.

  ! Local variables
  real :: t5((5*hi%iscb+1):(5*(hi%iecb+2)))  ! Temperatures along a line of subgrid locations [degC]
  real :: s5((5*hi%iscb+1):(5*(hi%iecb+2)))  ! Salinities along a line of subgrid locations [ppt]
  real :: t25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temperature variance along a line of subgrid locations [degC2]
  real :: ts5((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temp-salt covariance along a line of subgrid locations [degC ppt]
  real :: s25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS salinity variance along a line of subgrid locations [ppt2]
  real :: p5((5*hi%iscb+1):(5*(hi%iecb+2)))  ! Pressures along a line of subgrid locations, never
                                               ! rescaled from Pa [Pa]
  real :: r5((5*hi%iscb+1):(5*(hi%iecb+2)))  ! Densities anomalies along a line of subgrid
                                               ! locations [R ~> kg m-3] or [kg m-3]
  real :: t15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Temperatures at an array of subgrid locations [degC]
  real :: s15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Salinities at an array of subgrid locations [ppt]
  real :: t215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temperature variance along a line of subgrid locations [degC2]
  real :: ts15((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temp-salt covariance along a line of subgrid locations [degC ppt]
  real :: s215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS salinity variance along a line of subgrid locations [ppt2]
  real :: p15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Pressures at an array of subgrid locations [Pa]
  real :: r15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Densities at an array of subgrid locations
                                                 ! [R ~> kg m-3] or [kg m-3]
  real :: wt_t(5), wt_b(5)          ! Top and bottom weights [nondim]
  real :: rho_anom                  ! A density anomaly [R ~> kg m-3] or [kg m-3]
  real :: w_left, w_right           ! Left and right weights [nondim]
  real :: intz(5)    ! The gravitational acceleration times the integrals of density
                     ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]
  real, parameter :: c1_90 = 1.0/90.0  ! A rational constant [nondim]
  real :: gxrho      ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
  real :: i_rho      ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]
  real :: rho_scale  ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]
  real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]
  real :: dz(hi%iscb:hi%iecb+1)   ! Layer thicknesses at tracer points [Z ~> m]
  real :: dz_x(5,hi%iscb:hi%iecb) ! Layer thicknesses along an x-line of subrid locations [Z ~> m]
  real :: dz_y(5,hi%isc:hi%iec)   ! Layer thicknesses along a y-line of subrid locations [Z ~> m]
  real :: massweighttoggle          ! A non-dimensional toggle factor (0 or 1) [nondim]
  real :: ttl, tbl, ttr, tbr        ! Temperatures at the velocity cell corners [degC]
  real :: stl, sbl, str, sbr        ! Salinities at the velocity cell corners [ppt]
  real :: hwght                     ! A topographically limited thicknes weight [Z ~> m]
  real :: hl, hr                    ! Thicknesses to the left and right [Z ~> m]
  real :: idenom                    ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
  logical :: use_stanley_eos ! True is SGS variance fields exist in tv.
  logical :: use_vart, use_vars, use_covarts
  integer :: isq, ieq, jsq, jeq, i, j, m, n
  integer :: pos

  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB

  rho_scale = us%kg_m3_to_R
  gxrho = us%RL2_T2_to_Pa * g_e * rho_0
  rho_ref_mks = rho_ref * us%R_to_kg_m3
  i_rho = 1.0 / rho_0
  massweighttoggle = 0.
  if (present(usemasswghtinterp)) then
    if (usemasswghtinterp) massweighttoggle = 1.
  endif

  use_vart = associated(tv%varT)
  use_covarts = associated(tv%covarTS)
  use_vars = associated(tv%varS)
  use_stanley_eos = use_vart .or. use_covarts .or. use_vars
  t25(:) = 0.
  ts5(:) = 0.
  s25(:) = 0.
  t215(:) = 0.
  ts15(:) = 0.
  s215(:) = 0.

  do n = 1, 5
    wt_t(n) = 0.25 * real(5-n)
    wt_b(n) = 1.0 - wt_t(n)
  enddo

  ! 1. Compute vertical integrals
  do j=jsq,jeq+1
    do i = isq,ieq+1
      dz(i) = e(i,j,k) - e(i,j,k+1)
      do n=1,5
        p5(i*5+n) = -gxrho*(e(i,j,k) - 0.25*real(n-1)*dz(i))
        ! Salinity and temperature points are linearly interpolated
        s5(i*5+n) = wt_t(n) * s_t(i,j,k) + wt_b(n) * s_b(i,j,k)
        t5(i*5+n) = wt_t(n) * t_t(i,j,k) + wt_b(n) * t_b(i,j,k)
      enddo
      if (use_vart) t25(i*5+1:i*5+5) = tv%varT(i,j,k)
      if (use_covarts) ts5(i*5+1:i*5+5) = tv%covarTS(i,j,k)
      if (use_vars) s25(i*5+1:i*5+5) = tv%varS(i,j,k)
    enddo
    if (use_stanley_eos) then
      if (rho_scale /= 1.0) then
        call calculate_density(t5, s5, p5, t25, ts5, s25, r5, 1, (ieq-isq+2)*5, eos, &
                               rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t5, s5, p5, t25, ts5, s25, r5, 1, (ieq-isq+2)*5, eos, &
                               rho_ref=rho_ref_mks)
      endif
    else
      if (rho_scale /= 1.0) then
        call calculate_density(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho_ref=rho_ref_mks, &
                               scale=rho_scale)
      else
        call calculate_density(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho_ref=rho_ref_mks)
      endif
    endif

    do i=isq,ieq+1
    ! Use Boole's rule to estimate the pressure anomaly change.
      rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3))
      dpa(i,j) = g_e*dz(i)*rho_anom
      if (present(intz_dpa)) then
      ! Use a Boole's-rule-like fifth-order accurate estimate of
      ! the double integral of the pressure anomaly.
        intz_dpa(i,j) = 0.5*g_e*dz(i)**2 * &
                (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )
      endif
    enddo
  enddo ! end loops on j

  ! 2. Compute horizontal integrals in the x direction
  if (present(intx_dpa)) then ; do j=hi%jsc,hi%jec
    do i=isq,ieq
      ! Corner values of T and S
      ! hWght is the distance measure by which the cell is violation of
      ! hydrostatic consistency. For large hWght we bias the interpolation
      ! of T,S along the top and bottom integrals, almost like thickness
      ! weighting.
      ! Note: To work in terrain following coordinates we could offset
      ! this distance by the layer thickness to replicate other models.
      hwght = massweighttoggle * &
              max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k))
      if (hwght > 0.) then
        hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
        hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz_subroundoff
        hwght = hwght * ( (hl-hr)/(hl+hr) )**2
        idenom = 1./( hwght*(hr + hl) + hl*hr )
        ttl = ( (hwght*hr)*t_t(i+1,j,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
        ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i+1,j,k) ) * idenom
        tbl = ( (hwght*hr)*t_b(i+1,j,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
        tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i+1,j,k) ) * idenom
        stl = ( (hwght*hr)*s_t(i+1,j,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
        str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i+1,j,k) ) * idenom
        sbl = ( (hwght*hr)*s_b(i+1,j,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
        sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i+1,j,k) ) * idenom
      else
        ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i+1,j,k); tbr = t_b(i+1,j,k)
        stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i+1,j,k); sbr = s_b(i+1,j,k)
      endif

      do m=2,4
        w_left = wt_t(m) ; w_right = wt_b(m)
        dz_x(m,i) = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i+1,j,k) - e(i+1,j,k+1))

        ! Salinity and temperature points are linearly interpolated in
        ! the horizontal. The subscript (1) refers to the top value in
        ! the vertical profile while subscript (5) refers to the bottom
        ! value in the vertical profile.
        pos = i*15+(m-2)*5
        t15(pos+1) = w_left*ttl + w_right*ttr
        t15(pos+5) = w_left*tbl + w_right*tbr

        s15(pos+1) = w_left*stl + w_right*str
        s15(pos+5) = w_left*sbl + w_right*sbr

        p15(pos+1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i+1,j,k))

        ! Pressure
        do n=2,5
          p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)
        enddo

        ! Salinity and temperature (linear interpolation in the vertical)
        do n=2,4
          s15(pos+n) = wt_t(n) * s15(pos+1) + wt_b(n) * s15(pos+5)
          t15(pos+n) = wt_t(n) * t15(pos+1) + wt_b(n) * t15(pos+5)
        enddo
        if (use_vart) t215(pos+1:pos+5) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i+1,j,k)
        if (use_covarts) ts15(pos+1:pos+5) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i+1,j,k)
        if (use_vars) s215(pos+1:pos+5) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i+1,j,k)
      enddo
    enddo

    if (use_stanley_eos) then
      if (rho_scale /= 1.0) then
        call calculate_density(t15, s15, p15, t215, ts15, s215, r15, 1, 15*(ieq-isq+1), eos, &
                               rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t15, s15, p15, t215, ts15, s215, r15, 1, 15*(ieq-isq+1), eos, &
                               rho_ref=rho_ref_mks)
      endif
    else
      if (rho_scale /= 1.0) then
        call calculate_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho_ref=rho_ref_mks, &
                               scale=rho_scale)
      else
        call calculate_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho_ref=rho_ref_mks)
      endif
    endif

    do i=isq,ieq
      intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)

      ! Use Boole's rule to estimate the pressure anomaly change.
      do m = 2,4
        pos = i*15+(m-2)*5
        intz(m) = g_e*dz_x(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &
                          12.0*r15(pos+3)))
      enddo
      ! Use Boole's rule to integrate the bottom pressure anomaly values in x.
      intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                             12.0*intz(3))
    enddo
  enddo ; endif

  ! 3. Compute horizontal integrals in the y direction
  if (present(inty_dpa)) then ; do j=jsq,jeq
    do i=hi%isc,hi%iec
    ! Corner values of T and S
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, almost like thickness
    ! weighting.
    ! Note: To work in terrain following coordinates we could offset
    ! this distance by the layer thickness to replicate other models.
      hwght = massweighttoggle * &
              max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k))
      if (hwght > 0.) then
        hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
        hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz_subroundoff
        hwght = hwght * ( (hl-hr)/(hl+hr) )**2
        idenom = 1./( hwght*(hr + hl) + hl*hr )
        ttl = ( (hwght*hr)*t_t(i,j+1,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
        ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i,j+1,k) ) * idenom
        tbl = ( (hwght*hr)*t_b(i,j+1,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
        tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i,j+1,k) ) * idenom
        stl = ( (hwght*hr)*s_t(i,j+1,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
        str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i,j+1,k) ) * idenom
        sbl = ( (hwght*hr)*s_b(i,j+1,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
        sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i,j+1,k) ) * idenom
      else
        ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i,j+1,k); tbr = t_b(i,j+1,k)
        stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i,j+1,k); sbr = s_b(i,j+1,k)
      endif

      do m=2,4
        w_left = wt_t(m) ; w_right = wt_b(m)
        dz_y(m,i) = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i,j+1,k) - e(i,j+1,k+1))

        ! Salinity and temperature points are linearly interpolated in
        ! the horizontal. The subscript (1) refers to the top value in
        ! the vertical profile while subscript (5) refers to the bottom
        ! value in the vertical profile.
        pos = i*15+(m-2)*5
        t15(pos+1) = w_left*ttl + w_right*ttr
        t15(pos+5) = w_left*tbl + w_right*tbr

        s15(pos+1) = w_left*stl + w_right*str
        s15(pos+5) = w_left*sbl + w_right*sbr

        p15(pos+1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i,j+1,k))

        ! Pressure
        do n=2,5
          p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i)
        enddo

        ! Salinity and temperature (linear interpolation in the vertical)
        do n=2,4
          s15(pos+n) = wt_t(n) * s15(pos+1) + wt_b(n) * s15(pos+5)
          t15(pos+n) = wt_t(n) * t15(pos+1) + wt_b(n) * t15(pos+5)
        enddo
        if (use_vart) t215(pos+1:pos+5) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i,j+1,k)
        if (use_covarts) ts15(pos+1:pos+5) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i,j+1,k)
        if (use_vars) s215(pos+1:pos+5) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i,j+1,k)
      enddo
    enddo

    if (use_stanley_eos) then
      if (rho_scale /= 1.0) then
        call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
                               t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), &
                               r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, &
                               rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
                               t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), &
                               r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, rho_ref=rho_ref_mks)
      endif
    else
      if (rho_scale /= 1.0) then
        call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
                               r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, &
                               rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
                               r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, rho_ref=rho_ref_mks)
      endif
    endif

    do i=hi%isc,hi%iec
      intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)

      ! Use Boole's rule to estimate the pressure anomaly change.
      do m = 2,4
        pos = i*15+(m-2)*5
        intz(m) = g_e*dz_y(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &
                                         32.0*(r15(pos+2)+r15(pos+4)) + &
                                         12.0*r15(pos+3)))
      enddo
      ! Use Boole's rule to integrate the values.
      inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                             12.0*intz(3))
    enddo
  enddo ; endif

end subroutine int_density_dz_generic_plm


!> Compute pressure gradient force integrals for layer "k" and the case where T and S
!! are parabolic profiles
subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
                                      rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, &
                                      dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp)
  integer,              intent(in)  :: k   !< Layer index to calculate integrals for
  type(hor_index_type), intent(in)  :: hi  !< Ocean horizontal index structures for the input arrays
  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
  type(thermo_var_ptrs), intent(in) :: tv  !< Thermodynamic variables
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: t_t !< Potential temperature at the cell top [degC]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: t_b !< Potential temperature at the cell bottom [degC]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: s_t !< Salinity at the cell top [ppt]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                        intent(in)  :: s_b !< Salinity at the cell bottom [ppt]
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &
                        intent(in)  :: e   !< Height of interfaces [Z ~> m]
  real,                 intent(in)  :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is
                                           !! subtracted out to reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used to calculate
                                           !! the pressure (as p~=-z*rho_0*G_e) used in the equation of state.
  real,                 intent(in)  :: g_e !< The Earth's gravitational acceleration [m s-2]
  real,                 intent(in)  :: dz_subroundoff !< A minuscule thickness change [Z ~> m]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: bathyt !< The depth of the bathymetry [Z ~> m]
  type(eos_type),       pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in) :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [R L2 Z T-2 ~> Pa m]
  real, dimension(SZIB_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the x grid spacing [R L2 T-2 ~> Pa]
  real, dimension(SZI_(HI),SZJB_(HI)), &
              optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the y grid spacing [R L2 T-2 ~> Pa]
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting to
                                           !! interpolate T/S for top and bottom integrals.

! This subroutine calculates (by numerical quadrature) integrals of
! pressure anomalies across layers, which are required for calculating the
! finite-volume form pressure accelerations in a Boussinesq model.  The one
! potentially dodgy assumption here is that rho_0 is used both in the denominator
! of the accelerations, and in the pressure used to calculated density (the
! latter being -z*rho_0*G_e).  These two uses could be separated if need be.
!
! It is assumed that the salinity and temperature profiles are parabolic in the
! vertical. The top and bottom values within each layer are provided and
! a parabolic interpolation is used to compute intermediate values.

  ! Local variables
  real :: t5(5) ! Temperatures along a line of subgrid locations [degC]
  real :: s5(5) ! Salinities along a line of subgrid locations [ppt]
  real :: t25(5) ! SGS temperature variance along a line of subgrid locations [degC2]
  real :: ts5(5) ! SGS temperature-salinity covariance along a line of subgrid locations [degC ppt]
  real :: s25(5) ! SGS salinity variance along a line of subgrid locations [ppt2]
  real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa]
  real :: r5(5) ! Density anomalies from rho_ref at quadrature points [R ~> kg m-3] or [kg m-3]
  real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]
  real :: rho_anom ! The integrated density anomaly [R ~> kg m-3] or [kg m-3]
  real :: w_left, w_right  ! Left and right weights [nondim]
  real :: intz(5) ! The gravitational acceleration times the integrals of density
                  ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]
  real, parameter :: c1_90 = 1.0/90.0  ! Rational constants.
  real :: gxrho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
  real :: i_rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]
  real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]
  real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]
  real :: dz ! Layer thicknesses at tracer points [Z ~> m]
  real :: massweighttoggle ! A non-dimensional toggle factor (0 or 1) [nondim]
  real :: ttl, tbl, tml, ttr, tbr, tmr ! Temperatures at the velocity cell corners [degC]
  real :: stl, sbl, sml, str, sbr, smr ! Salinities at the velocity cell corners [ppt]
  real :: s6 ! PPM curvature coefficient for S [ppt]
  real :: t6 ! PPM curvature coefficient for T [degC]
  real :: t_top, t_mn, t_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T
  real :: s_top, s_mn, s_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S
  real :: hwght  ! A topographically limited thicknes weight [Z ~> m]
  real :: hl, hr ! Thicknesses to the left and right [Z ~> m]
  real :: idenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
  integer :: isq, ieq, jsq, jeq, i, j, m, n
  logical :: use_ppm ! If false, assume zero curvature in reconstruction, i.e. PLM
  logical :: use_stanley_eos ! True is SGS variance fields exist in tv.
  logical :: use_vart, use_vars, use_covarts

  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB

  rho_scale = us%kg_m3_to_R
  gxrho = us%RL2_T2_to_Pa * g_e * rho_0
  rho_ref_mks = rho_ref * us%R_to_kg_m3
  i_rho = 1.0 / rho_0
  massweighttoggle = 0.
  if (present(usemasswghtinterp)) then
    if (usemasswghtinterp) massweighttoggle = 1.
  endif

  ! In event PPM calculation is bypassed with use_PPM=False
  s6 = 0.
  t6 = 0.
  use_ppm = .true. ! This is a place-holder to allow later re-use of this function

  use_vart = associated(tv%varT)
  use_covarts = associated(tv%covarTS)
  use_vars = associated(tv%varS)
  use_stanley_eos = use_vart .or. use_covarts .or. use_vars
  t25(:) = 0.
  ts5(:) = 0.
  s25(:) = 0.

  do n = 1, 5
    wt_t(n) = 0.25 * real(5-n)
    wt_b(n) = 1.0 - wt_t(n)
  enddo

  ! 1. Compute vertical integrals
  do j=jsq,jeq+1 ; do i=isq,ieq+1
    if (use_ppm) then
      ! Curvature coefficient of the parabolas
      s6 = 3.0 * ( 2.0*tv%S(i,j,k) - ( s_t(i,j,k) + s_b(i,j,k) ) )
      t6 = 3.0 * ( 2.0*tv%T(i,j,k) - ( t_t(i,j,k) + t_b(i,j,k) ) )
    endif
    dz = e(i,j,k) - e(i,j,k+1)
    do n=1,5
      p5(n) = -gxrho*(e(i,j,k) - 0.25*real(n-1)*dz)
      ! Salinity and temperature points are reconstructed with PPM
      s5(n) = wt_t(n) * s_t(i,j,k) + wt_b(n) * ( s_b(i,j,k) + s6 * wt_t(n) )
      t5(n) = wt_t(n) * t_t(i,j,k) + wt_b(n) * ( t_b(i,j,k) + t6 * wt_t(n) )
    enddo
    if (use_stanley_eos) then
      if (use_vart) t25(:) = tv%varT(i,j,k)
      if (use_covarts) ts5(:) = tv%covarTS(i,j,k)
      if (use_vars) s25(:) = tv%varS(i,j,k)
      call calculate_density(t5, s5, p5, t25, ts5, s25, r5, &
                             1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
    else
      call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
    endif

    ! Use Boole's rule to estimate the pressure anomaly change.
    rho_anom = c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))
    dpa(i,j) = g_e*dz*rho_anom
    if (present(intz_dpa)) then
      ! Use a Boole's-rule-like fifth-order accurate estimate of
      ! the double integral of the pressure anomaly.
      intz_dpa(i,j) = 0.5*g_e*dz**2 * &
                      (rho_anom - c1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
    endif
  enddo ; enddo ! end loops on j and i

  ! 2. Compute horizontal integrals in the x direction
  if (present(intx_dpa)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
    ! Corner values of T and S
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, almost like thickness
    ! weighting.
    ! Note: To work in terrain following coordinates we could offset
    ! this distance by the layer thickness to replicate other models.
    hwght = massweighttoggle * &
            max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k))
    if (hwght > 0.) then
      hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
      hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz_subroundoff
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1./( hwght*(hr + hl) + hl*hr )
      ttl = ( (hwght*hr)*t_t(i+1,j,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
      tbl = ( (hwght*hr)*t_b(i+1,j,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
      tml = ( (hwght*hr)*tv%T(i+1,j,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom
      ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i+1,j,k) ) * idenom
      tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i+1,j,k) ) * idenom
      tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i+1,j,k) ) * idenom
      stl = ( (hwght*hr)*s_t(i+1,j,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
      sbl = ( (hwght*hr)*s_b(i+1,j,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
      sml = ( (hwght*hr)*tv%S(i+1,j,k) + (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom
      str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i+1,j,k) ) * idenom
      sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i+1,j,k) ) * idenom
      smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i+1,j,k) ) * idenom
    else
      ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i+1,j,k); tbr = t_b(i+1,j,k)
      tml = tv%T(i,j,k); tmr = tv%T(i+1,j,k)
      stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i+1,j,k); sbr = s_b(i+1,j,k)
      sml = tv%S(i,j,k); smr = tv%S(i+1,j,k)
    endif

    do m=2,4
      w_left = wt_t(m) ; w_right = wt_b(m)

      ! Salinity and temperature points are linearly interpolated in
      ! the horizontal. The subscript (1) refers to the top value in
      ! the vertical profile while subscript (5) refers to the bottom
      ! value in the vertical profile.
      t_top = w_left*ttl + w_right*ttr
      t_mn = w_left*tml + w_right*tmr
      t_bot = w_left*tbl + w_right*tbr

      s_top = w_left*stl + w_right*str
      s_mn = w_left*sml + w_right*smr
      s_bot = w_left*sbl + w_right*sbr

      ! Pressure
      dz = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i+1,j,k) - e(i+1,j,k+1))
      p5(1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i+1,j,k))
      do n=2,5
        p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo

      ! Parabolic reconstructions in the vertical for T and S
      if (use_ppm) then
        ! Coefficients of the parabolas
        s6 = 3.0 * ( 2.0*s_mn - ( s_top + s_bot ) )
        t6 = 3.0 * ( 2.0*t_mn - ( t_top + t_bot ) )
      endif
      do n=1,5
        s5(n) = wt_t(n) * s_top + wt_b(n) * ( s_bot + s6 * wt_t(n) )
        t5(n) = wt_t(n) * t_top + wt_b(n) * ( t_bot + t6 * wt_t(n) )
      enddo

      if (use_stanley_eos) then
        if (use_vart) t25(:) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i+1,j,k)
        if (use_covarts) ts5(:) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i+1,j,k)
        if (use_vars) s25(:) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i+1,j,k)
        call calculate_density(t5, s5, p5, t25, ts5, s25, r5, &
                               1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
      endif

      ! Use Boole's rule to estimate the pressure anomaly change.
      intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) )
    enddo ! m
    intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)

    ! Use Boole's rule to integrate the bottom pressure anomaly values in x.
    intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))

  enddo ; enddo ; endif

  ! 3. Compute horizontal integrals in the y direction
  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
    ! Corner values of T and S
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, almost like thickness
    ! weighting.
    ! Note: To work in terrain following coordinates we could offset
    ! this distance by the layer thickness to replicate other models.
    hwght = massweighttoggle * &
            max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k))
    if (hwght > 0.) then
      hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
      hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz_subroundoff
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1./( hwght*(hr + hl) + hl*hr )
      ttl = ( (hwght*hr)*t_t(i,j+1,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
      tbl = ( (hwght*hr)*t_b(i,j+1,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
      tml = ( (hwght*hr)*tv%T(i,j+1,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom
      ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i,j+1,k) ) * idenom
      tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i,j+1,k) ) * idenom
      tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i,j+1,k) ) * idenom
      stl = ( (hwght*hr)*s_t(i,j+1,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
      sbl = ( (hwght*hr)*s_b(i,j+1,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
      sml = ( (hwght*hr)*tv%S(i,j+1,k)+ (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom
      str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i,j+1,k) ) * idenom
      sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i,j+1,k) ) * idenom
      smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i,j+1,k) ) * idenom
    else
      ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i,j+1,k); tbr = t_b(i,j+1,k)
      tml = tv%T(i,j,k); tmr = tv%T(i,j+1,k)
      stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i,j+1,k); sbr = s_b(i,j+1,k)
      sml = tv%S(i,j,k); smr = tv%S(i,j+1,k)
    endif

    do m=2,4
      w_left = wt_t(m) ; w_right = wt_b(m)

      ! Salinity and temperature points are linearly interpolated in
      ! the horizontal. The subscript (1) refers to the top value in
      ! the vertical profile while subscript (5) refers to the bottom
      ! value in the vertical profile.
      t_top = w_left*ttl + w_right*ttr
      t_mn = w_left*tml + w_right*tmr
      t_bot = w_left*tbl + w_right*tbr

      s_top = w_left*stl + w_right*str
      s_mn = w_left*sml + w_right*smr
      s_bot = w_left*sbl + w_right*sbr

      ! Pressure
      dz = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i,j+1,k) - e(i,j+1,k+1))
      p5(1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i,j+1,k))
      do n=2,5
        p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo

      ! Parabolic reconstructions in the vertical for T and S
      if (use_ppm) then
        ! Coefficients of the parabolas
        s6 = 3.0 * ( 2.0*s_mn - ( s_top + s_bot ) )
        t6 = 3.0 * ( 2.0*t_mn - ( t_top + t_bot ) )
      endif
      do n=1,5
        s5(n) = wt_t(n) * s_top + wt_b(n) * ( s_bot + s6 * wt_t(n) )
        t5(n) = wt_t(n) * t_top + wt_b(n) * ( t_bot + t6 * wt_t(n) )
      enddo

      if (use_stanley_eos) then
        if (use_vart) t25(:) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i,j+1,k)
        if (use_covarts) ts5(:) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i,j+1,k)
        if (use_vars) s25(:) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i,j+1,k)
        call calculate_density(t5, s5, p5, t25, ts5, s25, r5, &
                               1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
      else
        call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
      endif

      ! Use Boole's rule to estimate the pressure anomaly change.
      intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) )
    enddo ! m
    intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)

    ! Use Boole's rule to integrate the bottom pressure anomaly values in y.
    inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))

  enddo ; enddo ; endif

end subroutine int_density_dz_generic_ppm

!> Calls the appropriate subroutine to calculate analytical and nearly-analytical
!! integrals in pressure across layers of geopotential anomalies, which are
!! required for calculating the finite-volume form pressure accelerations in a
!! non-Boussinesq model.  There are essentially no free assumptions, apart from the
!! use of Boole's rule to do the horizontal integrals, and from a truncation in the
!! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, &
                               dza, intp_dza, intx_dza, inty_dza, halo_size, &
                               bathyP, dP_tiny, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hi  !< The horizontal index structure
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: t   !< Potential temperature referenced to the surface [degC]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: s   !< Salinity [ppt]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: p_b !< Pressure at the bottom of the layer [R L2 T-2 ~> Pa] or [Pa]
  real,                 intent(in)  :: alpha_ref !< A mean specific volume that is subtracted out
                            !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
                            !! The calculation is mathematically identical with different values of
                            !! alpha_ref, but this reduces the effects of roundoff.
  type(eos_type),       pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in) :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(inout) :: dza !< The change in the geopotential anomaly across
                            !! the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the
                            !! geopotential anomaly relative to the anomaly at the bottom of the
                            !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]
  real, dimension(SZIB_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intx_dza !< The integral in x of the difference between the
                            !! geopotential anomaly at the top and bottom of the layer divided by
                            !! the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
  real, dimension(SZI_(HI),SZJB_(HI)), &
              optional, intent(inout) :: inty_dza !< The integral in y of the difference between the
                            !! geopotential anomaly at the top and bottom of the layer divided by
                            !! the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
  integer,    optional, intent(in)  :: halo_size !< The width of halo points on which to calculate dza.
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(in)  :: bathyp  !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
  real,       optional, intent(in)  :: dp_tiny !< A minuscule pressure change with
                            !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting
                            !! to interpolate T/S for top and bottom integrals.

  if (eos_quadrature(eos)) then
    call int_spec_vol_dp_generic_pcm(t, s, p_t, p_b, alpha_ref, hi, eos, us, &
                                     dza, intp_dza, intx_dza, inty_dza, halo_size, &
                                     bathyp, dp_tiny, usemasswghtinterp)
  else
    call analytic_int_specific_vol_dp(t, s, p_t, p_b, alpha_ref, hi, eos, &
                                      dza, intp_dza, intx_dza, inty_dza, halo_size, &
                                      bathyp, dp_tiny, usemasswghtinterp)
  endif

end subroutine int_specific_vol_dp


!>   This subroutine calculates integrals of specific volume anomalies in
!! pressure across layers, which are required for calculating the finite-volume
!! form pressure accelerations in a non-Boussinesq model.  There are essentially
!! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.
subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, dza, &
                                       intp_dza, intx_dza, inty_dza, halo_size, &
                                       bathyP, dP_neglect, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hi !< A horizontal index type structure.
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: t  !< Potential temperature of the layer [degC]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: s  !< Salinity of the layer [ppt]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]
  real,                 intent(in)  :: alpha_ref !< A mean specific volume that is subtracted out
                            !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
                            !! The calculation is mathematically identical with different values of
                            !! alpha_ref, but alpha_ref alters the effects of roundoff, and
                            !! answers do change.
  type(eos_type),       pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(inout) :: dza !< The change in the geopotential anomaly
                            !! across the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
                            !! the geopotential anomaly relative to the anomaly at the bottom of the
                            !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]
  real, dimension(SZIB_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intx_dza  !< The integral in x of the difference between
                            !! the geopotential anomaly at the top and bottom of the layer divided
                            !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
  real, dimension(SZI_(HI),SZJB_(HI)), &
              optional, intent(inout) :: inty_dza  !< The integral in y of the difference between
                            !! the geopotential anomaly at the top and bottom of the layer divided
                            !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
  integer,    optional, intent(in)  :: halo_size !< The width of halo points on which to calculate dza.
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(in)  :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
  real,       optional, intent(in)  :: dp_neglect !< A minuscule pressure change with
                                             !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting
                            !! to interpolate T/S for top and bottom integrals.

!   This subroutine calculates analytical and nearly-analytical integrals in
! pressure across layers of geopotential anomalies, which are required for
! calculating the finite-volume form pressure accelerations in a non-Boussinesq
! model.  There are essentially no free assumptions, apart from the use of
! Boole's rule to do the horizontal integrals, and from a truncation in the
! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.

  ! Local variables
  real :: t5(5)      ! Temperatures at five quadrature points [degC]
  real :: s5(5)      ! Salinities at five quadrature points [ppt]
  real :: p5(5)      ! Pressures at five quadrature points, scaled back to Pa if necessary [Pa]
  real :: a5(5)      ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]
  real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]
  real :: dp         ! The pressure change through a layer [R L2 T-2 ~> Pa]
  real :: hwght      ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
  real :: hl, hr     ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
  real :: alpha_ref_mks ! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1]
  real :: idenom     ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]
  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]
  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]
  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]
  real :: intp(5)    ! The integrals of specific volume with pressure at the
                     ! 5 sub-column locations [L2 T-2 ~> m2 s-2]
  real :: rl2_t2_to_pa  ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1]
  real :: sv_scale   ! A multiplicative factor by which to scale specific
                     ! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
  logical :: do_massweight ! Indicates whether to do mass weighting.
  real, parameter :: c1_90 = 1.0/90.0  ! A rational constant.
  integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, n, halo

  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
  halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
  ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
  if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh); endif
  if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh); endif

  sv_scale = us%R_to_kg_m3
  rl2_t2_to_pa = us%RL2_T2_to_Pa
  alpha_ref_mks = alpha_ref * us%kg_m3_to_R

  do_massweight = .false.
  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
    do_massweight = .true.
    if (.not.present(bathyp)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
        "bathyP must be present if useMassWghtInterp is present and true.")
    if (.not.present(dp_neglect)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
        "dP_neglect must be present if useMassWghtInterp is present and true.")
  endif ; endif

  do j=jsh,jeh ; do i=ish,ieh
    dp = p_b(i,j) - p_t(i,j)
    do n=1,5
      t5(n) = t(i,j) ; s5(n) = s(i,j)
      p5(n) = rl2_t2_to_pa * (p_b(i,j) - 0.25*real(n-1)*dp)
    enddo

    if (sv_scale /= 1.0) then
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
    else
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
    endif

    ! Use Boole's rule to estimate the interface height anomaly change.
    alpha_anom = c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3))
    dza(i,j) = dp*alpha_anom
    ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
    ! the interface height anomaly.
    if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
          (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
  enddo ; enddo

  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif

    intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr

      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness weighted.
      p5(1) = rl2_t2_to_pa * (wt_l*p_b(i,j) + wt_r*p_b(i+1,j))
      dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)

      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - rl2_t2_to_pa * 0.25*dp
      enddo
      if (sv_scale /= 1.0) then
        call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
      else
        call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
      endif

    ! Use Boole's rule to estimate the interface height anomaly change.
      intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
                                12.0*a5(3)))
    enddo
    ! Use Boole's rule to integrate the interface height anomaly values in x.
    intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
                           12.0*intp(3))
  enddo ; enddo ; endif

  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif

    intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr

      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness weighted.
      p5(1) = rl2_t2_to_pa * (wt_l*p_b(i,j) + wt_r*p_b(i,j+1))
      dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = rl2_t2_to_pa * (p5(n-1) - 0.25*dp)
      enddo
      if (sv_scale /= 1.0) then
        call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
      else
        call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
      endif

    ! Use Boole's rule to estimate the interface height anomaly change.
      intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
                                12.0*a5(3)))
    enddo
    ! Use Boole's rule to integrate the interface height anomaly values in y.
    inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
                           12.0*intp(3))
  enddo ; enddo ; endif

end subroutine int_spec_vol_dp_generic_pcm

!>   This subroutine calculates integrals of specific volume anomalies in
!! pressure across layers, which are required for calculating the finite-volume
!! form pressure accelerations in a non-Boussinesq model.  There are essentially
!! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.
subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, &
                             dP_neglect, bathyP, HI, EOS, US, dza, &
                             intp_dza, intx_dza, inty_dza, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hi !< A horizontal index type structure.
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: t_t  !< Potential temperature at the top of the layer [degC]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: t_b  !< Potential temperature at the bottom of the layer [degC]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: s_t  !< Salinity at the top the layer [ppt]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: s_b  !< Salinity at the bottom the layer [ppt]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]
  real,                 intent(in)  :: alpha_ref !< A mean specific volume that is subtracted out
                            !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
                            !! The calculation is mathematically identical with different values of
                            !! alpha_ref, but alpha_ref alters the effects of roundoff, and
                            !! answers do change.
  real,                 intent(in)  :: dp_neglect !<!< A miniscule pressure change with
                                             !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(in)  :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
  type(eos_type),       pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
  real, dimension(SZI_(HI),SZJ_(HI)), &
                        intent(inout) :: dza !< The change in the geopotential anomaly
                            !! across the layer [L2 T-2 ~> m2 s-2]
  real, dimension(SZI_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
                            !! the geopotential anomaly relative to the anomaly at the bottom of the
                            !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]
  real, dimension(SZIB_(HI),SZJ_(HI)), &
              optional, intent(inout) :: intx_dza  !< The integral in x of the difference between
                            !! the geopotential anomaly at the top and bottom of the layer divided
                            !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
  real, dimension(SZI_(HI),SZJB_(HI)), &
              optional, intent(inout) :: inty_dza  !< The integral in y of the difference between
                            !! the geopotential anomaly at the top and bottom of the layer divided
                            !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting
                            !! to interpolate T/S for top and bottom integrals.

!   This subroutine calculates analytical and nearly-analytical integrals in
! pressure across layers of geopotential anomalies, which are required for
! calculating the finite-volume form pressure accelerations in a non-Boussinesq
! model.  There are essentially no free assumptions, apart from the use of
! Boole's rule to do the horizontal integrals, and from a truncation in the
! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.

  real :: t5(5)      ! Temperatures at five quadrature points [degC]
  real :: s5(5)      ! Salinities at five quadrature points [ppt]
  real :: p5(5)      ! Pressures at five quadrature points, scaled back to Pa as necessary [Pa]
  real :: a5(5)      ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]
  real :: t15(15)    ! Temperatures at fifteen interior quadrature points [degC]
  real :: s15(15)    ! Salinities at fifteen interior quadrature points [ppt]
  real :: p15(15)    ! Pressures at fifteen quadrature points, scaled back to Pa as necessary [Pa]
  real :: a15(15)    ! Specific volumes at fifteen quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]
  real :: wt_t(5), wt_b(5) ! Weights of top and bottom values at quadrature points [nondim]
  real :: t_top, t_bot, s_top, s_bot, p_top, p_bot

  real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1]
  real :: dp         ! The pressure change through a layer [R L2 T-2 ~> Pa]
  real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa]
  real :: hwght      ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
  real :: hl, hr     ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
  real :: alpha_ref_mks ! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1]
  real :: idenom     ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]
  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]
  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]
  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]
  real :: intp(5)    ! The integrals of specific volume with pressure at the
                     ! 5 sub-column locations [L2 T-2 ~> m2 s-2]
  real :: rl2_t2_to_pa  ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1]
  real :: sv_scale   ! A multiplicative factor by which to scale specific
                     ! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
  real, parameter :: c1_90 = 1.0/90.0  ! A rational constant.
  logical :: do_massweight ! Indicates whether to do mass weighting.
  integer :: isq, ieq, jsq, jeq, i, j, m, n, pos

  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB

  do_massweight = .false.
  if (present(usemasswghtinterp)) do_massweight = usemasswghtinterp

  sv_scale = us%R_to_kg_m3
  rl2_t2_to_pa = us%RL2_T2_to_Pa
  alpha_ref_mks = alpha_ref * us%kg_m3_to_R

  do n = 1, 5 ! Note that these are reversed from int_density_dz.
    wt_t(n) = 0.25 * real(n-1)
    wt_b(n) = 1.0 - wt_t(n)
  enddo

  ! 1. Compute vertical integrals
  do j=jsq,jeq+1; do i=isq,ieq+1
    dp = p_b(i,j) - p_t(i,j)
    do n=1,5 ! T, S and p are linearly interpolated in the vertical.
      p5(n) = rl2_t2_to_pa * (wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j))
      s5(n) = wt_t(n) * s_t(i,j) + wt_b(n) * s_b(i,j)
      t5(n) = wt_t(n) * t_t(i,j) + wt_b(n) * t_b(i,j)
    enddo
    if (sv_scale /= 1.0) then
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
    else
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
    endif

    ! Use Boole's rule to estimate the interface height anomaly change.
    alpha_anom = c1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3))
    dza(i,j) = dp*alpha_anom
    ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
    ! the interface height anomaly.
    if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
          (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
  enddo ; enddo

  ! 2. Compute horizontal integrals in the x direction
  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, almost like thickness
    ! weighting. Note: To work in terrain following coordinates we could
    ! offset this distance by the layer thickness to replicate other models.
    hwght = 0.0
    if (do_massweight) &
      hwght =  max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif

    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr

      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness weighted.
      p_top = wt_l*p_t(i,j) + wt_r*p_t(i+1,j)
      p_bot = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
      t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i+1,j)
      t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i+1,j)
      s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i+1,j)
      s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i+1,j)
      dp_90(m) = c1_90*(p_bot - p_top)

      ! Salinity, temperature and pressure with linear interpolation in the vertical.
      pos = (m-2)*5
      do n=1,5
        p15(pos+n) = rl2_t2_to_pa * (wt_t(n) * p_top + wt_b(n) * p_bot)
        s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
        t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
      enddo
    enddo

    if (sv_scale /= 1.0) then
      call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks, scale=sv_scale)
    else
      call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks)
    endif

    intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
    do m=2,4
      ! Use Boole's rule to estimate the interface height anomaly change.
      ! The integrals at the ends of the segment are already known.
      pos = (m-2)*5
      intp(m) = dp_90(m)*((7.0*(a15(pos+1)+a15(pos+5)) + &
                          32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
    enddo
    ! Use Boole's rule to integrate the interface height anomaly values in x.
    intx_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
                           12.0*intp(3))
  enddo ; enddo ; endif

  ! 3. Compute horizontal integrals in the y direction
  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, like thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif

    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr

      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness weighted.
      p_top = wt_l*p_t(i,j) + wt_r*p_t(i,j+1)
      p_bot = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
      t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i,j+1)
      t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i,j+1)
      s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i,j+1)
      s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i,j+1)
      dp_90(m) = c1_90*(p_bot - p_top)

      ! Salinity, temperature and pressure with linear interpolation in the vertical.
      pos = (m-2)*5
      do n=1,5
        p15(pos+n) = rl2_t2_to_pa * (wt_t(n) * p_top + wt_b(n) * p_bot)
        s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
        t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
      enddo
    enddo

    if (sv_scale /= 1.0) then
      call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks, scale=sv_scale)
    else
      call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks)
    endif

    intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
    do m=2,4
      ! Use Boole's rule to estimate the interface height anomaly change.
      ! The integrals at the ends of the segment are already known.
      pos = (m-2)*5
      intp(m) = dp_90(m) * ((7.0*(a15(pos+1)+a15(pos+5)) + &
                            32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
    enddo
    ! Use Boole's rule to integrate the interface height anomaly values in x.
    inty_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
                           12.0*intp(3))
  enddo ; enddo ; endif

end subroutine int_spec_vol_dp_generic_plm


!> Find the depth at which the reconstructed pressure matches P_tgt
subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, &
                       rho_ref, G_e, EOS, US, P_b, z_out, z_tol)
  real,                  intent(in)  :: t_t !< Potential temperature at the cell top [degC]
  real,                  intent(in)  :: t_b !< Potential temperature at the cell bottom [degC]
  real,                  intent(in)  :: s_t !< Salinity at the cell top [ppt]
  real,                  intent(in)  :: s_b !< Salinity at the cell bottom [ppt]
  real,                  intent(in)  :: z_t !< Absolute height of top of cell [Z ~> m]   (Boussinesq ????)
  real,                  intent(in)  :: z_b !< Absolute height of bottom of cell [Z ~> m]
  real,                  intent(in)  :: p_t !< Anomalous pressure of top of cell, relative
                                            !! to g*rho_ref*z_t [R L2 T-2 ~> Pa]
  real,                  intent(in)  :: p_tgt !< Target pressure at height z_out, relative
                                            !! to g*rho_ref*z_out [R L2 T-2 ~> Pa]
  real,                  intent(in)  :: rho_ref !< Reference density with which calculation
                                            !! are anomalous to [R ~> kg m-3]
  real,                  intent(in)  :: g_e !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
  type(eos_type),        pointer     :: eos !< Equation of state structure
  type(unit_scale_type), intent(in)  :: us !< A dimensional unit scaling type
  real,                  intent(out) :: p_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa]
  real,                  intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m]
  real, optional,        intent(in)  :: z_tol !< The tolerance in finding z_out [Z ~> m]

  ! Local variables
  real :: dp    ! Pressure thickness of the layer [R L2 T-2 ~> Pa]
  real :: f_guess, f_l, f_r  ! Fractional positions [nondim]
  real :: gxrho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]
  real :: pa, pa_left, pa_right, pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz [R L2 T-2 ~> Pa]
  character(len=240) :: msg

  gxrho = g_e * rho_ref

  ! Anomalous pressure difference across whole cell
  dp = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, 1.0, eos)

  p_b = p_t + dp ! Anomalous pressure at bottom of cell

  if (p_tgt <= p_t ) then
    z_out = z_t
    return
  endif

  if (p_tgt >= p_b) then
    z_out = z_b
    return
  endif

  f_l = 0.
  pa_left = p_t - p_tgt ! Pa_left < 0
  f_r = 1.
  pa_right = p_b - p_tgt ! Pa_right > 0
  pa_tol = gxrho * 1.0e-5*us%m_to_Z
  if (present(z_tol)) pa_tol = gxrho * z_tol

  f_guess = f_l - pa_left / (pa_right - pa_left) * (f_r - f_l)
  pa = pa_right - pa_left ! To get into iterative loop
  do while ( abs(pa) > pa_tol )

    z_out = z_t + ( z_b - z_t ) * f_guess
    pa = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, f_guess, eos) - ( p_tgt - p_t )

    if (pa<pa_left) then
      write(msg,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
      call mom_error(fatal, 'find_depth_of_pressure_in_cell out of bounds negative: /n'//msg)
    elseif (pa<0.) then
      pa_left = pa
      f_l = f_guess
    elseif (pa>pa_right) then
      write(msg,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
      call mom_error(fatal, 'find_depth_of_pressure_in_cell out of bounds positive: /n'//msg)
    elseif (pa>0.) then
      pa_right = pa
      f_r = f_guess
    else ! Pa == 0
      return
    endif
    f_guess = f_l - pa_left / (pa_right - pa_left) * (f_r - f_l)

  enddo

end subroutine find_depth_of_pressure_in_cell


!> Returns change in anomalous pressure change from top to non-dimensional
!! position pos between z_t and z_b
real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS)
  real,           intent(in)  :: t_t !< Potential temperature at the cell top [degC]
  real,           intent(in)  :: t_b !< Potential temperature at the cell bottom [degC]
  real,           intent(in)  :: s_t !< Salinity at the cell top [ppt]
  real,           intent(in)  :: s_b !< Salinity at the cell bottom [ppt]
  real,           intent(in)  :: z_t !< The geometric height at the top of the layer [Z ~> m]
  real,           intent(in)  :: z_b !< The geometric height at the bottom of the layer [Z ~> m]
  real,           intent(in)  :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted out to
                                     !! reduce the magnitude of each of the integrals.
  real,           intent(in)  :: g_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
  real,           intent(in)  :: pos !< The fractional vertical position, 0 to 1 [nondim]
  type(eos_type), pointer     :: eos !< Equation of state structure
  real                        :: fract_dp_at_pos !< The change in pressure from the layer top to
                                     !! fractional position pos [R L2 T-2 ~> Pa]
  ! Local variables
  real, parameter :: c1_90 = 1.0/90.0  ! A rational constant [nondim]
  real :: dz                 ! Distance from the layer top [Z ~> m]
  real :: top_weight, bottom_weight ! Fractional weights at quadrature points [nondim]
  real :: rho_ave            ! Average density [R ~> kg m-3]
  real, dimension(5) :: t5   ! Temperatures at quadrature points [degC]
  real, dimension(5) :: s5   ! Salinities at quadrature points [ppt]
  real, dimension(5) :: p5   ! Pressures at quadrature points [R L2 T-2 ~> Pa]
  real, dimension(5) :: rho5 ! Densities at quadrature points [R ~> kg m-3]
  integer :: n

  do n=1,5
    ! Evaluate density at five quadrature points
    bottom_weight = 0.25*real(n-1) * pos
    top_weight = 1.0 - bottom_weight
    ! Salinity and temperature points are linearly interpolated
    s5(n) = top_weight * s_t + bottom_weight * s_b
    t5(n) = top_weight * t_t + bottom_weight * t_b
    p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( g_e * rho_ref )
  enddo
  call calculate_density(t5, s5, p5, rho5, eos)
  rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref

  ! Use Boole's rule to estimate the average density
  rho_ave = c1_90*(7.0*(rho5(1)+rho5(5)) + 32.0*(rho5(2)+rho5(4)) + 12.0*rho5(3))

  dz = ( z_t - z_b ) * pos
  frac_dp_at_pos = g_e * dz * rho_ave
end function frac_dp_at_pos

end module mom_density_integrals

!> \namespace mom_density_integrals
!!