MOM_EOS.F90

!> Provides subroutines for quantities specific to the equation of state
module mom_eos
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_eos_linear, only : calculate_density_linear, calculate_spec_vol_linear
use mom_eos_linear, only : calculate_density_derivs_linear
use mom_eos_linear, only : calculate_specvol_derivs_linear, int_density_dz_linear
use mom_eos_linear, only : calculate_density_second_derivs_linear
use mom_eos_linear, only : calculate_compress_linear, int_spec_vol_dp_linear
use mom_eos_wright, only : calculate_density_wright, calculate_spec_vol_wright
use mom_eos_wright, only : calculate_density_derivs_wright
use mom_eos_wright, only : calculate_specvol_derivs_wright, int_density_dz_wright
use mom_eos_wright, only : calculate_compress_wright, int_spec_vol_dp_wright
use mom_eos_wright, only : calculate_density_second_derivs_wright
use mom_eos_unesco, only : calculate_density_unesco, calculate_spec_vol_unesco
use mom_eos_unesco, only : calculate_density_derivs_unesco, calculate_density_unesco
use mom_eos_unesco, only : calculate_compress_unesco
use mom_eos_nemo,   only : calculate_density_nemo
use mom_eos_nemo,   only : calculate_density_derivs_nemo, calculate_density_nemo
use mom_eos_nemo,   only : calculate_compress_nemo
use mom_eos_teos10, only : calculate_density_teos10, calculate_spec_vol_teos10
use mom_eos_teos10, only : calculate_density_derivs_teos10
use mom_eos_teos10, only : calculate_specvol_derivs_teos10
use mom_eos_teos10, only : calculate_density_second_derivs_teos10
use mom_eos_teos10, only : calculate_compress_teos10
use mom_eos_teos10, only : gsw_sp_from_sr, gsw_pt_from_ct
use mom_tfreeze,    only : calculate_tfreeze_linear, calculate_tfreeze_millero
use mom_tfreeze,    only : calculate_tfreeze_teos10
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_string_functions, only : uppercase
use mom_hor_index, only : hor_index_type
 
implicit none ; private
 
#include <MOM_memory.h>
 
public calculate_compress, calculate_density, query_compressible
public calculate_density_derivs, calculate_specific_vol_derivs
public calculate_density_second_derivs
public eos_init, eos_manual_init, eos_end, eos_allocate
public eos_use_linear, calculate_spec_vol
public int_density_dz, int_specific_vol_dp
public int_density_dz_generic_plm, int_density_dz_generic_ppm
public int_spec_vol_dp_generic_plm !, int_spec_vol_dz_generic_ppm
public int_density_dz_generic, int_spec_vol_dp_generic
public find_depth_of_pressure_in_cell
public calculate_tfreeze
public convert_temp_salt_for_teos10
public gsw_sp_from_sr, gsw_pt_from_ct
public extract_member_eos
 
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
 
!> Calculates density of sea water from T, S and P
interface calculate_density
  module procedure calculate_density_scalar, calculate_density_array
end interface calculate_density
 
!> Calculates specific volume of sea water from T, S and P
interface calculate_spec_vol
  module procedure calculate_spec_vol_scalar, calculate_spec_vol_array
end interface calculate_spec_vol
 
!> Calculate the derivatives of density with temperature and salinity from T, S, and P
interface calculate_density_derivs
  module procedure calculate_density_derivs_scalar, calculate_density_derivs_array
end interface calculate_density_derivs
 
!> Calculates the second derivatives of density with various combinations of temperature,
!! salinity, and pressure from T, S and P
interface calculate_density_second_derivs
  module procedure calculate_density_second_derivs_scalar, calculate_density_second_derivs_array
end interface calculate_density_second_derivs
 
!> Calculates the freezing point of sea water from T, S and P
interface calculate_tfreeze
  module procedure calculate_tfreeze_scalar, calculate_tfreeze_array
end interface calculate_tfreeze
 
!> A control structure for the equation of state
type, public :: eos_type ; private
  integer :: form_of_eos = 0 !< The equation of state to use.
  integer :: form_of_tfreeze = 0 !< The expression for the potential temperature
                             !! of the freezing point.
  logical :: eos_quadrature  !< If true, always use the generic (quadrature)
                             !! code for the integrals of density.
  logical :: compressible = .true. !< If true, in situ density is a function of pressure.
! The following parameters are used with the linear equation of state only.
  real :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
  real :: drho_dt   !< The partial derivative of density with temperature [kg m-3 degC-1]
  real :: drho_ds   !< The partial derivative of density with salinity [kg m-3 ppt-1].
! The following parameters are use with the linear expression for the freezing
! point only.
  real :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC].
  real :: dtfr_ds   !< The derivative of freezing point with salinity [degC ppt-1].
  real :: dtfr_dp   !< The derivative of freezing point with pressure [degC Pa-1].
 
!  logical :: test_EOS = .true. ! If true, test the equation of state
end type eos_type
 
! The named integers that might be stored in eqn_of_state_type%form_of_EOS.
integer, parameter, public :: eos_linear = 1 !< A named integer specifying an equation of state
integer, parameter, public :: eos_unesco = 2 !< A named integer specifying an equation of state
integer, parameter, public :: eos_wright = 3 !< A named integer specifying an equation of state
integer, parameter, public :: eos_teos10 = 4 !< A named integer specifying an equation of state
integer, parameter, public :: eos_nemo   = 5 !< A named integer specifying an equation of state
 
character*(10), parameter :: eos_linear_string = "LINEAR" !< A string for specifying the equation of state
character*(10), parameter :: eos_unesco_string = "UNESCO" !< A string for specifying the equation of state
character*(10), parameter :: eos_wright_string = "WRIGHT" !< A string for specifying the equation of state
character*(10), parameter :: eos_teos10_string = "TEOS10" !< A string for specifying the equation of state
character*(10), parameter :: eos_nemo_string   = "NEMO"   !< A string for specifying the equation of state
character*(10), parameter :: eos_default = eos_wright_string !< The default equation of state
 
integer, parameter :: tfreeze_linear = 1  !< A named integer specifying a freezing point expression
integer, parameter :: tfreeze_millero = 2 !< A named integer specifying a freezing point expression
integer, parameter :: tfreeze_teos10 = 3  !< A named integer specifying a freezing point expression
character*(10), parameter :: tfreeze_linear_string = "LINEAR" !< A string for specifying the freezing point expression
character*(10), parameter :: tfreeze_millero_string = "MILLERO_78" !< A string for specifying
                                                              !! freezing point expression
character*(10), parameter :: tfreeze_teos10_string = "TEOS10" !< A string for specifying the freezing point expression
character*(10), parameter :: tfreeze_default = tfreeze_linear_string !< The default freezing point expression
 
contains
 
!> Calls the appropriate subroutine to calculate density of sea water for scalar inputs.
!! If rho_ref is present, the anomaly with respect to rho_ref is returned.
subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref)
  real,           intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3]
  type(eos_type), pointer     :: EOS !< Equation of state structure
  real, optional, intent(in)  :: rho_ref  !< A reference density [kg m-3].
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_linear(t, s, pressure, rho, &
                                      eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
    case (eos_unesco)
      call calculate_density_unesco(t, s, pressure, rho, rho_ref)
    case (eos_wright)
      call calculate_density_wright(t, s, pressure, rho, rho_ref)
    case (eos_teos10)
      call calculate_density_teos10(t, s, pressure, rho, rho_ref)
    case (eos_nemo)
      call calculate_density_nemo(t, s, pressure, rho, rho_ref)
    case default
      call mom_error(fatal, &
           "calculate_density_scalar: EOS is not valid.")
  end select
 
end subroutine calculate_density_scalar
 
!> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs.
!! If rho_ref is present, the anomaly with respect to rho_ref is returned.
subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_ref)
  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3]
  integer,            intent(in)  :: start !< Start index for computation
  integer,            intent(in)  :: npts !< Number of point to compute
  type(eos_type),     pointer     :: EOS !< Equation of state structure
  real,     optional, intent(in)  :: rho_ref  !< A reference density [kg m-3].
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_array called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_linear(t, s, pressure, rho, start, npts, &
                                      eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
    case (eos_unesco)
      call calculate_density_unesco(t, s, pressure, rho, start, npts, rho_ref)
    case (eos_wright)
      call calculate_density_wright(t, s, pressure, rho, start, npts, rho_ref)
    case (eos_teos10)
      call calculate_density_teos10(t, s, pressure, rho, start, npts, rho_ref)
    case (eos_nemo)
      call calculate_density_nemo  (t, s, pressure, rho, start, npts, rho_ref)
    case default
      call mom_error(fatal, &
           "calculate_density_array: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_density_array
 
!> Calls the appropriate subroutine to calculate specific volume of sea water
!! for scalar inputs.
subroutine calculate_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref)
  real,           intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: specvol  !< specific volume (in-situ if pressure is local) [m3 kg-1]
  type(eos_type), pointer     :: EOS      !< Equation of state structure
  real, optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
 
  real :: rho
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_spec_vol_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_spec_vol_linear(t, s, pressure, specvol, &
               eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
    case (eos_unesco)
      call calculate_spec_vol_unesco(t, s, pressure, specvol, spv_ref)
    case (eos_wright)
      call calculate_spec_vol_wright(t, s, pressure, specvol, spv_ref)
    case (eos_teos10)
      call calculate_spec_vol_teos10(t, s, pressure, specvol, spv_ref)
    case (eos_nemo)
      call calculate_density_nemo(t, s, pressure, rho)
      if (present(spv_ref)) then
        specvol = 1.0 / rho - spv_ref
      else
        specvol = 1.0 / rho
      endif
    case default
      call mom_error(fatal, &
           "calculate_spec_vol_scalar: EOS is not valid.")
  end select
 
end subroutine calculate_spec_vol_scalar
 
 
!> Calls the appropriate subroutine to calculate the specific volume of sea water
!! for 1-D array inputs.
subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref)
  real, dimension(:), intent(in)  :: T        !< potential temperature relative to the surface
                                              !! [degC].
  real, dimension(:), intent(in)  :: S        !< salinity [ppt].
  real, dimension(:), intent(in)  :: pressure !< pressure [Pa].
  real, dimension(:), intent(out) :: specvol  !< in situ specific volume [kg m-3].
  integer,            intent(in)  :: start    !< the starting point in the arrays.
  integer,            intent(in)  :: npts     !< the number of values to calculate.
  type(eos_type),     pointer     :: EOS      !< Equation of state structure
  real,     optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
 
  real, dimension(size(specvol)) :: rho
 
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_spec_vol_array called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_spec_vol_linear(t, s, pressure, specvol, start, npts, &
               eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
    case (eos_unesco)
      call calculate_spec_vol_unesco(t, s, pressure, specvol, start, npts, spv_ref)
    case (eos_wright)
      call calculate_spec_vol_wright(t, s, pressure, specvol, start, npts, spv_ref)
    case (eos_teos10)
      call calculate_spec_vol_teos10(t, s, pressure, specvol, start, npts, spv_ref)
    case (eos_nemo)
      call calculate_density_nemo  (t, s, pressure, rho, start, npts)
      if (present(spv_ref)) then
        specvol(:) = 1.0 / rho(:) - spv_ref
      else
        specvol(:) = 1.0 / rho(:)
      endif
    case default
      call mom_error(fatal, &
           "calculate_spec_vol_array: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_spec_vol_array
 
 
!> Calls the appropriate subroutine to calculate the freezing point for scalar inputs.
subroutine calculate_tfreeze_scalar(S, pressure, T_fr, EOS)
  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: T_fr !< Freezing point potential temperature referenced
                                      !! to the surface [degC]
  type(eos_type), pointer     :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_TFreeze)
    case (tfreeze_linear)
      call calculate_tfreeze_linear(s, pressure, t_fr, eos%TFr_S0_P0, &
                                    eos%dTFr_dS, eos%dTFr_dp)
    case (tfreeze_millero)
      call calculate_tfreeze_millero(s, pressure, t_fr)
    case (tfreeze_teos10)
      call calculate_tfreeze_teos10(s, pressure, t_fr)
    case default
      call mom_error(fatal, &
           "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
  end select
 
end subroutine calculate_tfreeze_scalar
 
!> Calls the appropriate subroutine to calculate the freezing point for a 1-D array.
subroutine calculate_tfreeze_array(S, pressure, T_fr, start, npts, EOS)
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: T_fr !< Freezing point potential temperature referenced
                                          !! to the surface [degC]
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(eos_type),     pointer     :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_TFreeze)
    case (tfreeze_linear)
      call calculate_tfreeze_linear(s, pressure, t_fr, start, npts, &
                                    eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
    case (tfreeze_millero)
      call calculate_tfreeze_millero(s, pressure, t_fr, start, npts)
    case (tfreeze_teos10)
      call calculate_tfreeze_teos10(s, pressure, t_fr, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
  end select
 
end subroutine calculate_tfreeze_array
 
!> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.
subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS)
  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: drho_dT !< The partial derivative of density with potential
                                             !! temperature [kg m-3 degC-1].
  real, dimension(:), intent(out) :: drho_dS !< The partial derivative of density with salinity,
                                             !! in [kg m-3 ppt-1].
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(eos_type),     pointer     :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_derivs_linear(t, s, pressure, drho_dt, drho_ds, eos%Rho_T0_S0, &
                                           eos%dRho_dT, eos%dRho_dS, start, npts)
    case (eos_unesco)
      call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
    case (eos_wright)
      call calculate_density_derivs_wright(t, s, pressure, drho_dt, drho_ds, start, npts)
    case (eos_teos10)
      call calculate_density_derivs_teos10(t, s, pressure, drho_dt, drho_ds, start, npts)
    case (eos_nemo)
      call calculate_density_derivs_nemo(t, s, pressure, drho_dt, drho_ds, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs_array: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_density_derivs_array
 
!> Calls the appropriate subroutines to calculate density derivatives by promoting a scalar
!! to a one-element array
subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS)
  real,           intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: drho_dT !< The partial derivative of density with potential
                                         !! temperature [kg m-3 degC-1].
  real,           intent(out) :: drho_dS !< The partial derivative of density with salinity,
                                         !! in [kg m-3 ppt-1].
  type(eos_type), pointer     :: EOS !< Equation of state structure
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_derivs_linear(t, s, pressure, drho_dt, drho_ds, &
                                                  eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
    case (eos_wright)
      call calculate_density_derivs_wright(t, s, pressure, drho_dt, drho_ds)
    case (eos_teos10)
      call calculate_density_derivs_teos10(t, s, pressure, drho_dt, drho_ds)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs_scalar: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_density_derivs_scalar
 
!> Calls the appropriate subroutine to calculate density second derivatives for 1-D array inputs.
subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, &
                                                 drho_dS_dP, drho_dT_dP, start, npts, EOS)
  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect
                                                !! to S [kg m-3 ppt-2]
  real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with respcct
                                                !! to T [kg m-3 ppt-1 degC-1]
  real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect
                                                !! to T [kg m-3 degC-2]
  real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect
                                                !! to pressure [kg m-3 ppt-1 Pa-1]
  real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect
                                                !! to pressure [kg m-3 degC-1 Pa-1]
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(eos_type),     pointer     :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_second_derivs_linear(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
    case (eos_wright)
      call calculate_density_second_derivs_wright(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
    case (eos_teos10)
      call calculate_density_second_derivs_teos10(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_density_second_derivs_array
 
!> Calls the appropriate subroutine to calculate density second derivatives for scalar nputs.
subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, &
                                                  drho_dS_dP, drho_dT_dP, EOS)
  real, intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, intent(in)  :: S !< Salinity [ppt]
  real, intent(in)  :: pressure !< Pressure [Pa]
  real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect
                                  !! to S [kg m-3 ppt-2]
  real, intent(out) :: drho_dS_dT !< Partial derivative of beta with respcct
                                  !! to T [kg m-3 ppt-1 degC-1]
  real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect
                                  !! to T [kg m-3 degC-2]
  real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect
                                  !! to pressure [kg m-3 ppt-1 Pa-1]
  real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect
                                  !! to pressure [kg m-3 degC-1 Pa-1]
  type(eos_type),     pointer     :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_second_derivs_linear(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp)
    case (eos_wright)
      call calculate_density_second_derivs_wright(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp)
    case (eos_teos10)
      call calculate_density_second_derivs_teos10(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_density_second_derivs_scalar
 
!> Calls the appropriate subroutine to calculate specific volume derivatives for an array.
subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS)
  real, dimension(:), intent(in)  :: t !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: s !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: dsv_dt !< The partial derivative of specific volume with potential
                                            !! temperature [m3 kg-1 degC-1].
  real, dimension(:), intent(out) :: dsv_ds !< The partial derivative of specific volume with salinity
                                            !! [m3 kg-1 ppt-1].
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(eos_type),     pointer     :: eos !< Equation of state structure
  ! Local variables
  real, dimension(size(T)) :: drho_dt, drho_ds, rho
  integer :: j
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_specvol_derivs_linear(t, s, pressure, dsv_dt, dsv_ds, start, &
                                           npts, eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
    case (eos_unesco)
      call calculate_density_unesco(t, s, pressure, rho, start, npts)
      call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
      do j=start,start+npts-1
        dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
        dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
      enddo
    case (eos_wright)
      call calculate_specvol_derivs_wright(t, s, pressure, dsv_dt, dsv_ds, start, npts)
    case (eos_teos10)
      call calculate_specvol_derivs_teos10(t, s, pressure, dsv_dt, dsv_ds, start, npts)
    case (eos_nemo)
      call calculate_density_nemo(t, s, pressure, rho, start, npts)
      call calculate_density_derivs_nemo(t, s, pressure, drho_dt, drho_ds, start, npts)
      do j=start,start+npts-1
        dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
        dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
      enddo
    case default
      call mom_error(fatal, &
           "calculate_density_derivs: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_specific_vol_derivs
 
!> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs.
subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS)
  real, dimension(:), intent(in)  :: t !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: s !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: rho !< In situ density [kg m-3].
  real, dimension(:), intent(out) :: drho_dp !< The partial derivative of density with pressure
                                     !! (also the inverse of the square of sound speed) in s2 m-2.
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(eos_type),     pointer     :: eos !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_compress called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_compress_linear(t, s, pressure, rho, drho_dp, start, npts, &
                                     eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
    case (eos_unesco)
      call calculate_compress_unesco(t, s, pressure, rho, drho_dp, start, npts)
    case (eos_wright)
      call calculate_compress_wright(t, s, pressure, rho, drho_dp, start, npts)
    case (eos_teos10)
      call calculate_compress_teos10(t, s, pressure, rho, drho_dp, start, npts)
    case (eos_nemo)
      call calculate_compress_nemo(t, s, pressure, rho, drho_dp, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_compress: EOS%form_of_EOS is not valid.")
  end select
 
end subroutine calculate_compress
 
!> Calls the appropriate subroutine to alculate 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 Bode's rule to do the horizontal integrals, and from a truncation in the
!! series for log(1-eps/1+eps) that assumes that |eps| <  .
subroutine 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)
  type(hor_index_type), intent(in)  :: hi  !< The horizontal index structure
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: t   !< Potential temperature referenced to the surface [degC]
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: s   !< Salinity [ppt]
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_t !< Pressure at the top of the layer [Pa].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_b !< Pressure at the bottom of the layer [Pa].
  real,                 intent(in)  :: alpha_ref !< A mean specific volume that is subtracted out
                            !! to reduce the magnitude of each of the integrals, 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
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(out) :: dza !< The change in the geopotential anomaly across
                            !! the layer [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(out) :: intp_dza !< The integral in pressure through the layer of the
                            !! geopotential anomaly relative to the anomaly at the bottom of the
                            !! layer [Pa m2 s-2].
  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
              optional, intent(out) :: 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 [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
              optional, intent(out) :: 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 [m2 s-2].
  integer,    optional, intent(in)  :: halo_size !< The width of halo points on which to calculate dza.
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(in)  :: bathyp  !< The pressure at the bathymetry [Pa]
  real,       optional, intent(in)  :: dp_tiny !< A miniscule pressure change with
                                               !! the same units as p_t (Pa?)
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting
                            !! to interpolate T/S for top and bottom integrals.
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "int_specific_vol_dp called with an unassociated EOS_type EOS.")
 
  if (eos%EOS_quadrature) then
    call int_spec_vol_dp_generic(t, s, p_t, p_b, alpha_ref, hi, eos, &
                                 dza, intp_dza, intx_dza, inty_dza, halo_size, &
                                 bathyp, dp_tiny, usemasswghtinterp)
  else ; select case (eos%form_of_EOS)
    case (eos_linear)
      call int_spec_vol_dp_linear(t, s, p_t, p_b, alpha_ref, hi, eos%Rho_T0_S0, &
                                  eos%dRho_dT, eos%dRho_dS, dza, intp_dza, &
                                  intx_dza, inty_dza, halo_size, &
                                  bathyp, dp_tiny, usemasswghtinterp)
    case (eos_wright)
      call int_spec_vol_dp_wright(t, s, p_t, p_b, alpha_ref, hi, dza, &
                                  intp_dza, intx_dza, inty_dza, halo_size, &
                                  bathyp, dp_tiny, usemasswghtinterp)
    case default
      call int_spec_vol_dp_generic(t, s, p_t, p_b, alpha_ref, hi, eos, &
                                   dza, intp_dza, intx_dza, inty_dza, halo_size, &
                                   bathyp, dp_tiny, usemasswghtinterp)
  end select ; endif
 
end subroutine int_specific_vol_dp
 
!> This subroutine calculates analytical and nearly-analytical 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(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, &
                          dpa, intz_dpa, intx_dpa, inty_dpa, &
                          bathyT, dz_neglect, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hii !< Ocean horizontal index structures for the input arrays
  type(hor_index_type), intent(in)  :: hio !< Ocean horizontal index structures for the output arrays
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: t   !< Potential temperature referenced to the surface [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: s   !< Salinity [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< Height at the top of the layer in depth units [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is subtracted out to
                                           !! reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [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 [m2 Z-1 s-2 ~> m s-2].
  type(eos_type),       pointer     :: eos !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [Pa Z ~> Pa m].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: 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 [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: 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 [Pa].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
              optional, intent(in)  :: bathyt !< The depth of the bathymetry [Z ~> m].
  real,       optional, intent(in)  :: dz_neglect !< A miniscule thickness change [Z ~> m].
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting to
                                           !! interpolate T/S for top and bottom integrals.
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "int_density_dz called with an unassociated EOS_type EOS.")
 
  if (eos%EOS_quadrature) then
    call int_density_dz_generic(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
                                eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
                                bathyt, dz_neglect, usemasswghtinterp)
  else ; select case (eos%form_of_EOS)
    case (eos_linear)
      call int_density_dz_linear(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio,  &
                                 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, &
                                 dpa, intz_dpa, intx_dpa, inty_dpa, &
                                 bathyt, dz_neglect, usemasswghtinterp)
    case (eos_wright)
      call int_density_dz_wright(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio,  &
                                 dpa, intz_dpa, intx_dpa, inty_dpa, &
                                 bathyt, dz_neglect, usemasswghtinterp)
    case default
      call int_density_dz_generic(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio,  &
                                  eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
                                  bathyt, dz_neglect, usemasswghtinterp)
  end select ; endif
 
end subroutine int_density_dz
 
!> Returns true if the equation of state is compressible (i.e. has pressure dependence)
logical function query_compressible(EOS)
  type(eos_type), pointer :: eos !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "query_compressible called with an unassociated EOS_type EOS.")
 
  query_compressible = eos%compressible
end function query_compressible
 
!> Initializes EOS_type by allocating and reading parameters
subroutine eos_init(param_file, EOS)
  type(param_file_type), intent(in) :: param_file !< Parameter file structure
  type(eos_type),        pointer    :: eos !< Equation of state structure
  ! Local variables
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_EOS" ! This module's name.
  character(len=40)  :: tmpstr
 
  if (.not.associated(eos)) call eos_allocate(eos)
 
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
 
  call get_param(param_file, mdl, "EQN_OF_STATE", tmpstr, &
                 "EQN_OF_STATE determines which ocean equation of state "//&
                 "should be used.  Currently, the valid choices are "//&
                 '"LINEAR", "UNESCO", "WRIGHT", "NEMO" and "TEOS10". '//&
                 "This is only used if USE_EOS is true.", default=eos_default)
  select case (uppercase(tmpstr))
    case (eos_linear_string)
      eos%form_of_EOS = eos_linear
    case (eos_unesco_string)
      eos%form_of_EOS = eos_unesco
    case (eos_wright_string)
      eos%form_of_EOS = eos_wright
    case (eos_teos10_string)
      eos%form_of_EOS = eos_teos10
    case (eos_nemo_string)
      eos%form_of_EOS = eos_nemo
    case default
      call mom_error(fatal, "interpret_eos_selection: EQN_OF_STATE "//&
                              trim(tmpstr) // "in input file is invalid.")
  end select
  call mom_mesg('interpret_eos_selection: equation of state set to "' // &
                trim(tmpstr)//'"', 5)
 
  if (eos%form_of_EOS == eos_linear) then
    eos%Compressible = .false.
    call get_param(param_file, mdl, "RHO_T0_S0", eos%Rho_T0_S0, &
                 "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
                 "this is the density at T=0, S=0.", units="kg m-3", &
                 default=1000.0)
    call get_param(param_file, mdl, "DRHO_DT", eos%dRho_dT, &
                 "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
                 "this is the partial derivative of density with "//&
                 "temperature.", units="kg m-3 K-1", default=-0.2)
    call get_param(param_file, mdl, "DRHO_DS", eos%dRho_dS, &
                 "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
                 "this is the partial derivative of density with "//&
                 "salinity.", units="kg m-3 PSU-1", default=0.8)
  endif
 
  call get_param(param_file, mdl, "EOS_QUADRATURE", eos%EOS_quadrature, &
                 "If true, always use the generic (quadrature) code "//&
                 "code for the integrals of density.", default=.false.)
 
  call get_param(param_file, mdl, "TFREEZE_FORM", tmpstr, &
                 "TFREEZE_FORM determines which expression should be "//&
                 "used for the freezing point.  Currently, the valid "//&
                 'choices are "LINEAR", "MILLERO_78", "TEOS10"', &
                 default=tfreeze_default)
  select case (uppercase(tmpstr))
    case (tfreeze_linear_string)
      eos%form_of_TFreeze = tfreeze_linear
    case (tfreeze_millero_string)
      eos%form_of_TFreeze = tfreeze_millero
    case (tfreeze_teos10_string)
      eos%form_of_TFreeze = tfreeze_teos10
    case default
      call mom_error(fatal, "interpret_eos_selection:  TFREEZE_FORM "//&
                              trim(tmpstr) // "in input file is invalid.")
  end select
 
  if (eos%form_of_TFreeze == tfreeze_linear) then
    call get_param(param_file, mdl, "TFREEZE_S0_P0",eos%TFr_S0_P0, &
                 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
                 "this is the freezing potential temperature at "//&
                 "S=0, P=0.", units="deg C", default=0.0)
    call get_param(param_file, mdl, "DTFREEZE_DS",eos%dTFr_dS, &
                 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
                 "this is the derivative of the freezing potential "//&
                 "temperature with salinity.", &
                 units="deg C PSU-1", default=-0.054)
    call get_param(param_file, mdl, "DTFREEZE_DP",eos%dTFr_dP, &
                 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
                 "this is the derivative of the freezing potential "//&
                 "temperature with pressure.", &
                 units="deg C Pa-1", default=0.0)
  endif
 
  if ((eos%form_of_EOS == eos_teos10 .OR. eos%form_of_EOS == eos_nemo) .AND. &
      eos%form_of_TFreeze /= tfreeze_teos10) then
      call mom_error(fatal, "interpret_eos_selection:  EOS_TEOS10 or EOS_NEMO \n" //&
      "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 .")
  endif
 
 
end subroutine eos_init
 
!> Manually initialized an EOS type (intended for unit testing of routines which need a specific EOS)
subroutine eos_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, &
                           Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
  type(eos_type),    pointer    :: eos !< Equation of state structure
  integer, optional, intent(in) :: form_of_eos !< A coded integer indicating the equation of state to use.
  integer, optional, intent(in) :: form_of_tfreeze !< A coded integer indicating the expression for
                                       !! the potential temperature of the freezing point.
  logical, optional, intent(in) :: eos_quadrature !< If true, always use the generic (quadrature)
                                       !! code for the integrals of density.
  logical, optional, intent(in) :: compressible  !< If true, in situ density is a function of pressure.
  real   , optional, intent(in) :: rho_t0_s0 !< Density at T=0 degC and S=0 ppt [kg m-3]
  real   , optional, intent(in) :: drho_dt   !< Partial derivative of density with temperature
                                             !! in [kg m-3 degC-1]
  real   , optional, intent(in) :: drho_ds   !< Partial derivative of density with salinity
                                             !! in [kg m-3 ppt-1]
  real   , optional, intent(in) :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC].
  real   , optional, intent(in) :: dtfr_ds   !< The derivative of freezing point with salinity
                                             !! in [degC ppt-1].
  real   , optional, intent(in) :: dtfr_dp   !< The derivative of freezing point with pressure
                                             !! in [degC Pa-1].
 
  if (present(form_of_eos    ))  eos%form_of_EOS     = form_of_eos
  if (present(form_of_tfreeze))  eos%form_of_TFreeze = form_of_tfreeze
  if (present(eos_quadrature ))  eos%EOS_quadrature  = eos_quadrature
  if (present(compressible   ))  eos%Compressible    = compressible
  if (present(rho_t0_s0      ))  eos%Rho_T0_S0       = rho_t0_s0
  if (present(drho_dt        ))  eos%drho_dT         = drho_dt
  if (present(drho_ds        ))  eos%dRho_dS         = drho_ds
  if (present(tfr_s0_p0      ))  eos%TFr_S0_P0       = tfr_s0_p0
  if (present(dtfr_ds        ))  eos%dTFr_dS         = dtfr_ds
  if (present(dtfr_dp        ))  eos%dTFr_dp         = dtfr_dp
 
end subroutine eos_manual_init
 
!> Allocates EOS_type
subroutine eos_allocate(EOS)
  type(eos_type), pointer :: eos !< Equation of state structure
 
  if (.not.associated(eos)) allocate(eos)
end subroutine eos_allocate
 
!> Deallocates EOS_type
subroutine eos_end(EOS)
  type(eos_type), pointer :: eos !< Equation of state structure
 
  if (associated(eos)) deallocate(eos)
end subroutine eos_end
 
!> Set equation of state structure (EOS) to linear with given coefficients
!!
!! \note This routine is primarily for testing and allows a local copy of the
!! EOS_type (EOS argument) to be set to use the linear equation of state
!! independent from the rest of the model.
subroutine eos_use_linear(Rho_T0_S0, dRho_dT, dRho_dS, EOS, use_quadrature)
  real,              intent(in) :: rho_t0_s0 !< Density at T=0 degC and S=0 ppt [kg m-3]
  real,              intent(in) :: drho_dt   !< Partial derivative of density with temperature [kg m-3 degC-1]
  real,              intent(in) :: drho_ds   !< Partial derivative of density with salinity [kg m-3 ppt-1]
  logical, optional, intent(in) :: use_quadrature !< If true, always use the generic (quadrature)
                                             !! code for the integrals of density.
  type(eos_type),    pointer    :: eos       !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "MOM_EOS.F90: EOS_use_linear() called with an unassociated EOS_type EOS.")
 
  eos%form_of_EOS = eos_linear
  eos%Compressible = .false.
  eos%Rho_T0_S0 = rho_t0_s0
  eos%dRho_dT = drho_dt
  eos%dRho_dS = drho_ds
  eos%EOS_quadrature = .false.
  if (present(use_quadrature)) eos%EOS_quadrature = use_quadrature
 
end subroutine eos_use_linear
 
!>   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.
subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
                                  EOS, dpa, intz_dpa, intx_dpa, inty_dpa, &
                                  bathyT, dz_neglect, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hii !< Horizontal index type for input variables.
  type(hor_index_type), intent(in)  :: hio !< Horizontal index type for output variables.
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: t  !< Potential temperature of the layer [degC].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: s  !< Salinity of the layer [ppt].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< Height at the top of the layer in depth units [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is
                                          !! subtracted out to reduce the magnitude
                                          !! of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [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 [m2 Z-1 s-2 ~> m s-2].
  type(eos_type),       pointer     :: eos !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly
                                          !! across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the
                                          !! layer of the pressure anomaly relative to the
                                          !! anomaly at the top of the layer [Pa Z ~> Pa m].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: 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 [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: 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 [Pa].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
              optional, intent(in)  :: bathyt !< The depth of the bathymetry [Z ~> m].
  real,       optional, intent(in)  :: dz_neglect !< A miniscule thickness change [Z ~> m].
  logical,    optional, intent(in)  :: usemasswghtinterp !< If true, uses mass weighting to
                                          !! interpolate T/S for top and bottom integrals.
  real :: t5(5), s5(5), p5(5), r5(5)
  real :: rho_anom   ! The depth averaged density anomaly [kg m-3].
  real :: w_left, w_right
  real, parameter :: c1_90 = 1.0/90.0  ! Rational constants.
  real :: gxrho, i_rho
  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 [Pa].
  logical :: do_massweight ! Indicates whether to do mass weighting.
  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n, ioff, joff
 
  ioff = hio%idg_offset - hii%idg_offset
  joff = hio%jdg_offset - hii%jdg_offset
 
  ! These array bounds work for the indexing convention of the input arrays, but
  ! on the computational domain defined for the output arrays.
  isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
  jsq = hio%JscB + joff ; jeq = hio%JecB + joff
  is = hio%isc + ioff ; ie = hio%iec + ioff
  js = hio%jsc + joff ; je = hio%jec + joff
 
  gxrho = g_e * rho_0
  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
    call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
 
    ! Use Bode'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-ioff,j-joff) = g_e*dz*rho_anom
    ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of
    ! the pressure anomaly.
    if (present(intz_dpa)) intz_dpa(i-ioff,j-joff) = 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-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
    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
      call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
 
    ! Use Bode'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 Bode's rule to integrate the bottom pressure anomaly values in x.
    intx_dpa(i-ioff,j-joff) = 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-ioff,j-joff) ; intz(5) = dpa(i-ioff,j-joff+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
      call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
 
    ! Use Bode'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 Bode's rule to integrate the values.
    inty_dpa(i-ioff,j-joff) = 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
 
 
! ==========================================================================
!> Compute pressure gradient force integrals by quadrature for the case where
!! T and S are linear profiles.
subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, &
                                       rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, &
                                       intz_dpa, intx_dpa, inty_dpa, &
                                       useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hii !< Ocean horizontal index structures for the input arrays
  type(hor_index_type), intent(in)  :: hio !< Ocean horizontal index structures for the output arrays
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: t_t !< Potential temperatue at the cell top [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: t_b !< Potential temperatue at the cell bottom [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: s_t !< Salinity at the cell top [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: s_b !< Salinity at the cell bottom [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< The geometric height at the top of the layer,
                                           !! in depth units [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< The geometric height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is subtracted out to
                                           !! reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [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 [m2 Z-1 s-2 ~> m s-2].
  real,                 intent(in)  :: dz_subroundoff !< A miniscule thickness change [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: bathyt !< The depth of the bathymetry [Z ~> m].
  type(eos_type),       pointer     :: eos !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [Pa Z].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: 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 [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: 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 [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 assumtion 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*hio%iscb+1):(5*(hio%iecb+2)))  ! Temperatures along a line of subgrid locations [degC].
  real :: s5((5*hio%iscb+1):(5*(hio%iecb+2)))  ! Salinities along a line of subgrid locations [ppt].
  real :: p5((5*hio%iscb+1):(5*(hio%iecb+2)))  ! Pressures along a line of subgrid locations [Pa].
  real :: r5((5*hio%iscb+1):(5*(hio%iecb+2)))  ! Densities along a line of subgrid locations [kg m-3].
  real :: t15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Temperatures at an array of subgrid locations [degC].
  real :: s15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Salinities at an array of subgrid locations [ppt].
  real :: p15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Pressures at an array of subgrid locations [Pa].
  real :: r15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Densities at an array of subgrid locations [kg m-3].
  real :: wt_t(5), wt_b(5)          ! Top and bottom weights [nondim].
  real :: rho_anom                  ! A density anomaly [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 [Pa].
  real, parameter :: c1_90 = 1.0/90.0  ! A rational constant [nondim].
  real :: gxrho                     ! Gravitational acceleration times density [kg m-1 Z-1 s-2 ~> kg m-2 s-2].
  real :: i_rho                     ! The inverse of the reference density [m3 kg-1].
  real :: dz(hio%iscb:hio%iecb+1)   ! Layer thicknesses at tracer points [Z ~> m].
  real :: dz_x(5,hio%iscb:hio%iecb) ! Layer thicknesses along an x-line of subrid locations [Z ~> m].
  real :: dz_y(5,hio%isc:hio%iec)   ! Layer thicknesses along a y-line of subrid locations [Z ~> m].
  real :: weight_t, weight_b        ! Nondimensional weights of the top and bottom.
  real :: massweighttoggle          ! A nondimensional toggle factor (0 or 1).
  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].
  integer :: isq, ieq, jsq, jeq, i, j, m, n
  integer :: iin, jin, ioff, joff
  integer :: pos
 
  ioff = hio%idg_offset - hii%idg_offset
  joff = hio%jdg_offset - hii%jdg_offset
 
  isq = hio%IscB ; ieq = hio%IecB ; jsq = hio%JscB ; jeq = hio%JecB
 
  gxrho = g_e * rho_0
  i_rho = 1.0 / rho_0
  massweighttoggle = 0.
  if (present(usemasswghtinterp)) then
    if (usemasswghtinterp) massweighttoggle = 1.
  endif
 
  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
    jin = j+joff
    do i = isq,ieq+1 ; iin = i+ioff
      dz(i) = z_t(iin,jin) - z_b(iin,jin)
      do n=1,5
        p5(i*5+n) = -gxrho*(z_t(iin,jin) - 0.25*real(n-1)*dz(i))
        ! Salinity and temperature points are linearly interpolated
        s5(i*5+n) = wt_t(n) * s_t(iin,jin) + wt_b(n) * s_b(iin,jin)
        t5(i*5+n) = wt_t(n) * t_t(iin,jin) + wt_b(n) * t_b(iin,jin)
      enddo
    enddo
    call calculate_density_array(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho_ref )
 
    do i=isq,ieq+1 ; iin = i+ioff
    ! Use Bode'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 Bode'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=hio%jsc,hio%jec ; jin = j+joff
    do i=isq,ieq ; iin = i+ioff
      ! 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(iin,jin)-z_t(iin+1,jin), -bathyt(iin+1,jin)-z_t(iin,jin))
      if (hwght > 0.) then
        hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
        hr = (z_t(iin+1,jin) - z_b(iin+1,jin)) + dz_subroundoff
        hwght = hwght * ( (hl-hr)/(hl+hr) )**2
        idenom = 1./( hwght*(hr + hl) + hl*hr )
        ttl = ( (hwght*hr)*t_t(iin+1,jin) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
        ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin+1,jin) ) * idenom
        tbl = ( (hwght*hr)*t_b(iin+1,jin) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
        tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin+1,jin) ) * idenom
        stl = ( (hwght*hr)*s_t(iin+1,jin) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
        str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin+1,jin) ) * idenom
        sbl = ( (hwght*hr)*s_b(iin+1,jin) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
        sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin+1,jin) ) * idenom
      else
        ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin+1,jin); tbr = t_b(iin+1,jin)
        stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin+1,jin); sbr = s_b(iin+1,jin)
      endif
 
      do m=2,4
        w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
        dz_x(m,i) = w_left*(z_t(iin,jin) - z_b(iin,jin)) + w_right*(z_t(iin+1,jin) - z_b(iin+1,jin))
 
        ! 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*z_t(iin,jin) + w_right*z_t(iin+1,jin))
 
        ! 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
          weight_t = 0.25 * real(5-n)
          weight_b = 1.0 - weight_t
          s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
          t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
        enddo
      enddo
    enddo
 
    call calculate_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho_ref)
 
    do i=isq,ieq ; iin = i+ioff
      intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
 
      ! Use Bode'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 Bode'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 ; jin = j+joff
    do i=hio%isc,hio%iec ; iin = i+ioff
    ! 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)-z_t(iin,jin+1), -bathyt(i,j+1)-z_t(iin,jin))
      if (hwght > 0.) then
        hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
        hr = (z_t(iin,jin+1) - z_b(iin,jin+1)) + dz_subroundoff
        hwght = hwght * ( (hl-hr)/(hl+hr) )**2
        idenom = 1./( hwght*(hr + hl) + hl*hr )
        ttl = ( (hwght*hr)*t_t(iin,jin+1) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
        ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin,jin+1) ) * idenom
        tbl = ( (hwght*hr)*t_b(iin,jin+1) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
        tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin,jin+1) ) * idenom
        stl = ( (hwght*hr)*s_t(iin,jin+1) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
        str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin,jin+1) ) * idenom
        sbl = ( (hwght*hr)*s_b(iin,jin+1) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
        sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin,jin+1) ) * idenom
      else
        ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin,jin+1); tbr = t_b(iin,jin+1)
        stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin,jin+1); sbr = s_b(iin,jin+1)
      endif
 
      do m=2,4
        w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
        dz_y(m,i) = w_left*(z_t(iin,jin) - z_b(iin,jin)) + w_right*(z_t(iin,jin+1) - z_b(iin,jin+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*z_t(iin,jin) + w_right*z_t(iin,jin+1))
 
        ! 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
          weight_t = 0.25 * real(5-n)
          weight_b = 1.0 - weight_t
          s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
          t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
        enddo
      enddo
    enddo
 
    call calculate_density_array(t15(15*hio%isc+1:), s15(15*hio%isc+1:), p15(15*hio%isc+1:), &
                                 r15(15*hio%isc+1:), 1, 15*(hio%iec-hio%isc+1), eos, rho_ref)
    do i=hio%isc,hio%iec ; iin = i+ioff
      intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
 
      ! Use Bode'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 Bode'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
! ==========================================================================
! Above is the routine where only the S and T profiles are modified
! The real topography is still used
! ==========================================================================
 
!> 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, P_b, z_out, z_tol)
  real,           intent(in)  :: t_t !< Potential temperatue at the cell top [degC]
  real,           intent(in)  :: t_b !< Potential temperatue 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 [Pa]
  real,           intent(in)  :: p_tgt !< Target pressure at height z_out, relative to g*rho_ref*z_out [Pa]
  real,           intent(in)  :: rho_ref !< Reference density with which calculation are anomalous to
  real,           intent(in)  :: g_e !< Gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
  type(eos_type), pointer     :: eos !< Equation of state structure
  real,           intent(out) :: p_b !< Pressure at the bottom of the cell [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 :: top_weight, bottom_weight, rho_anom, w_left, w_right, gxrho, dz, dp, f_guess, f_l, f_r
  real :: pa, pa_left, pa_right, pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz
 
  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.e-5 ! 1e-5 has dimensions of m, but should be converted to the units of 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(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
      stop 'Blurgh! Too negative'
    elseif (pa<0.) then
      pa_left = pa
      f_l = f_guess
    elseif (pa>pa_right) then
      write(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
      stop 'Blurgh! Too positive'
    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 temperatue at the cell top [degC]
  real,           intent(in)  :: t_b !< Potential temperatue 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 [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 [m s-2]
  real,           intent(in)  :: pos !< The fractional vertical position, 0 to 1 [nondim].
  type(eos_type), pointer     :: eos !< Equation of state structure
  ! Local variables
  real, parameter :: c1_90 = 1.0/90.0  ! Rational constants.
  real :: dz, top_weight, bottom_weight, rho_ave
  real, dimension(5) :: t5, s5, p5, rho5
  integer :: n
 
  do n=1,5
    ! Evalute 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_array(t5, s5, p5, rho5, 1, 5, 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
 
 
! ==========================================================================
!> Compute pressure gradient force integrals for the case where T and S
!! are parabolic profiles
subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, &
                                       z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
                                       EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
 
  type(hor_index_type), intent(in)  :: hii !< Ocean horizontal index structures for the input arrays
  type(hor_index_type), intent(in)  :: hio !< Ocean horizontal index structures for the output arrays
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: t   !< Potential temperature referenced to the surface [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: t_t !< Potential temperatue at the cell top [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: t_b !< Potential temperatue at the cell bottom [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: s   !< Salinity [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: s_t !< Salinity at the cell top [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: s_b !< Salinity at the cell bottom [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< Height at the top of the layer [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is subtracted out to
                                           !! reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [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]
  type(eos_type),       pointer     :: eos !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [Pa Z ~> Pa m].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: 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 [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: 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 [Pa].
 
! 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 assumtion 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), s5(5), p5(5), r5(5)
  real :: rho_anom
  real :: w_left, w_right, intz(5)
  real, parameter :: c1_90 = 1.0/90.0  ! Rational constants.
  real :: gxrho, i_rho
  real :: dz
  real :: weight_t, weight_b
  real :: s0, s1, s2                   ! parabola coefficients for S [ppt]
  real :: t0, t1, t2                   ! parabola coefficients for T [degC]
  real :: xi                           ! normalized coordinate
  real :: t_top, t_mid, t_bot
  real :: s_top, s_mid, s_bot
  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n, ioff, joff
  real, dimension(4) :: x, y
  real, dimension(9) :: s_node, t_node, p_node, r_node
 
 
  call mom_error(fatal, &
    "int_density_dz_generic_ppm: the implementation is not done yet, contact developer")
 
  ioff = hio%idg_offset - hii%idg_offset
  joff = hio%jdg_offset - hii%jdg_offset
 
  ! These array bounds work for the indexing convention of the input arrays, but
  ! on the computational domain defined for the output arrays.
  isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
  jsq = hio%JscB + joff ; jeq = hio%JecB + joff
  is = hio%isc + ioff ; ie = hio%iec + ioff
  js = hio%jsc + joff ; je = hio%jec + joff
 
  gxrho = g_e * rho_0
  i_rho = 1.0 / rho_0
 
  ! =============================
  ! 1. Compute vertical integrals
  ! =============================
  do j=jsq,jeq+1 ; do i=isq,ieq+1
    dz = z_t(i,j) - z_b(i,j)
 
    ! Coefficients of the parabola for S
    s0 = s_t(i,j)
    s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
    s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0*s(i,j) )
 
    ! Coefficients of the parabola for T
    t0 = t_t(i,j)
    t1 = 6.0 * t(i,j) - 4.0 * t_t(i,j) - 2.0 * t_b(i,j)
    t2 = 3.0 * ( t_t(i,j) + t_b(i,j) - 2.0*t(i,j) )
 
    do n=1,5
      p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
 
      ! Parabolic reconstruction for T and S
      xi = 0.25 * ( n - 1 )
      s5(n) = s0 + s1 * xi + s2 * xi**2
      t5(n) = t0 + t1 * xi + t2 * xi**2
    enddo
 
    call calculate_density(t5, s5, p5, r5, 1, 5, eos)
 
    ! Use Bode'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)) - &
    !       rho_ref
 
    rho_anom = 1000.0 + s(i,j) - rho_ref
    dpa(i-ioff,j-joff) = g_e*dz*rho_anom
 
    ! Use a Bode's-rule-like fifth-order accurate estimate of
    ! the double integral of the pressure anomaly.
    !r5 = r5 - rho_ref
    !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))) )
 
    intz_dpa(i-ioff,j-joff) = 0.5 * g_e * dz**2 * ( 1000.0 - rho_ref + s0 + s1/3.0 + &
                                    s2/6.0 )
  enddo ; enddo ! end loops on j and i
 
  ! ==================================================
  ! 2. Compute horizontal integrals in the x direction
  ! ==================================================
  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
    intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
    do m=2,4
      w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
      dz = w_left*(z_t(i,j) - z_b(i,j)) + w_right*(z_t(i+1,j) - z_b(i+1,j))
 
      ! 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*t_t(i,j) + w_right*t_t(i+1,j)
      t_mid = w_left*t(i,j)   + w_right*t(i+1,j)
      t_bot = w_left*t_b(i,j) + w_right*t_b(i+1,j)
 
      s_top = w_left*s_t(i,j) + w_right*s_t(i+1,j)
      s_mid = w_left*s(i,j)   + w_right*s(i+1,j)
      s_bot = w_left*s_b(i,j) + w_right*s_b(i+1,j)
 
      p5(1) = -gxrho*(w_left*z_t(i,j) + w_right*z_t(i+1,j))
 
      ! Pressure
      do n=2,5
        p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo
 
      ! Coefficients of the parabola for S
      s0 = s_top
      s1 = 6.0 * s_mid - 4.0 * s_top - 2.0 * s_bot
      s2 = 3.0 * ( s_top + s_bot - 2.0*s_mid )
 
      ! Coefficients of the parabola for T
      t0 = t_top
      t1 = 6.0 * t_mid - 4.0 * t_top - 2.0 * t_bot
      t2 = 3.0 * ( t_top + t_bot - 2.0*t_mid )
 
      do n=1,5
        ! Parabolic reconstruction for T and S
        xi = 0.25 * ( n - 1 )
        s5(n) = s0 + s1 * xi + s2 * xi**2
        t5(n) = t0 + t1 * xi + t2 * xi**2
      enddo
 
      call calculate_density(t5, s5, p5, r5, 1, 5, eos)
 
    ! Use Bode'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)) - rho_ref)
    enddo
    intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                           12.0*intz(3))
 
    ! Use Gauss quadrature rule to compute integral
 
    ! The following coordinates define the quadrilateral on which the integral
    ! is computed
    x(1) = 1.0
    x(2) = 0.0
    x(3) = 0.0
    x(4) = 1.0
    y(1) = z_t(i+1,j)
    y(2) = z_t(i,j)
    y(3) = z_b(i,j)
    y(4) = z_b(i+1,j)
 
    t_node = 0.0
    p_node = 0.0
 
    ! Nodal values for S
 
    ! Parabolic reconstruction on the left
    s0 = s_t(i,j)
    s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
    s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0 * s(i,j) )
    s_node(2) = s0
    s_node(6) = s0 + 0.5 * s1 + 0.25 * s2
    s_node(3) = s0 + s1 + s2
 
    ! Parabolic reconstruction on the left
    s0 = s_t(i+1,j)
    s1 = 6.0 * s(i+1,j) - 4.0 * s_t(i+1,j) - 2.0 * s_b(i+1,j)
    s2 = 3.0 * ( s_t(i+1,j) + s_b(i+1,j) - 2.0 * s(i+1,j) )
    s_node(1) = s0
    s_node(8) = s0 + 0.5 * s1 + 0.25 * s2
    s_node(4) = s0 + s1 + s2
 
    s_node(5) = 0.5 * ( s_node(2) + s_node(1) )
    s_node(9) = 0.5 * ( s_node(6) + s_node(8) )
    s_node(7) = 0.5 * ( s_node(3) + s_node(4) )
 
    call calculate_density( t_node, s_node, p_node, r_node, 1, 9, eos )
    r_node = r_node - rho_ref
 
    call compute_integral_quadratic( x, y, r_node, intx_dpa(i-ioff,j-joff) )
 
    intx_dpa(i-ioff,j-joff) = intx_dpa(i-ioff,j-joff) * g_e
 
  enddo ; enddo ; endif
 
  ! ==================================================
  ! 3. Compute horizontal integrals in the y direction
  ! ==================================================
  if (present(inty_dpa)) then
    call mom_error(warning, "int_density_dz_generic_ppm still needs to be written for inty_dpa!")
    do j=jsq,jeq ; do i=is,ie
 
      inty_dpa(i-ioff,j-joff) = 0.0
 
    enddo ; enddo
  endif
 
end subroutine int_density_dz_generic_ppm
 
 
 
! =============================================================================
!> Compute the integral of the quadratic function
subroutine compute_integral_quadratic( x, y, f, integral )
  real, dimension(4), intent(in)  :: x  !< The x-position of the corners
  real, dimension(4), intent(in)  :: y  !< The y-position of the corners
  real, dimension(9), intent(in)  :: f  !< The function at the quadrature points
  real,               intent(out) :: integral !< The returned integral
 
  ! Local variables
  integer               :: i, k
  real, dimension(9)    :: weight, xi, eta          ! integration points
  real                  :: f_k
  real                  :: dxdxi, dxdeta
  real                  :: dydxi, dydeta
  real, dimension(4)    :: phiiso, dphiisodxi, dphiisodeta
  real, dimension(9)    :: phi, dphidxi, dphideta
  real                  :: jacobian_k
  real                  :: t
 
  ! Quadrature rule (4 points)
  !weight(:) = 1.0
  !xi(1) = - sqrt(3.0) / 3.0
  !xi(2) = sqrt(3.0) / 3.0
  !xi(3) = sqrt(3.0) / 3.0
  !xi(4) = - sqrt(3.0) / 3.0
  !eta(1) = - sqrt(3.0) / 3.0
  !eta(2) = - sqrt(3.0) / 3.0
  !eta(3) = sqrt(3.0) / 3.0
  !eta(4) = sqrt(3.0) / 3.0
 
  ! Quadrature rule (9 points)
  t = sqrt(3.0/5.0)
  weight(1) = 25.0/81.0 ; xi(1) = -t ; eta(1) = t
  weight(2) = 40.0/81.0 ; xi(2) = .0 ; eta(2) = t
  weight(3) = 25.0/81.0 ; xi(3) =  t ; eta(3) = t
  weight(4) = 40.0/81.0 ; xi(4) = -t ; eta(4) = .0
  weight(5) = 64.0/81.0 ; xi(5) = .0 ; eta(5) = .0
  weight(6) = 40.0/81.0 ; xi(6) =  t ; eta(6) = .0
  weight(7) = 25.0/81.0 ; xi(7) = -t ; eta(7) = -t
  weight(8) = 40.0/81.0 ; xi(8) = .0 ; eta(8) = -t
  weight(9) = 25.0/81.0 ; xi(9) =  t ; eta(9) = -t
 
  integral = 0.0
 
  ! Integration loop
  do k = 1,9
 
    ! Evaluate shape functions and gradients for isomorphism
    call evaluate_shape_bilinear( xi(k), eta(k), phiiso, &
                                  dphiisodxi, dphiisodeta )
 
    ! Determine gradient of global coordinate at integration point
    dxdxi  = 0.0
    dxdeta = 0.0
    dydxi  = 0.0
    dydeta = 0.0
 
    do i = 1,4
      dxdxi  = dxdxi  + x(i) * dphiisodxi(i)
      dxdeta = dxdeta + x(i) * dphiisodeta(i)
      dydxi  = dydxi  + y(i) * dphiisodxi(i)
      dydeta = dydeta + y(i) * dphiisodeta(i)
    enddo
 
    ! Evaluate Jacobian at integration point
    jacobian_k = dxdxi*dydeta - dydxi*dxdeta
 
    ! Evaluate shape functions for interpolation
    call evaluate_shape_quadratic( xi(k), eta(k), phi, dphidxi, dphideta )
 
    ! Evaluate function at integration point
    f_k = 0.0
    do i = 1,9
      f_k = f_k + f(i) * phi(i)
    enddo
 
    integral = integral + weight(k) * f_k * jacobian_k
 
  enddo ! end integration loop
 
end subroutine compute_integral_quadratic
 
 
! =============================================================================
!> Evaluation of the four bilinear shape fn and their gradients at (xi,eta)
subroutine evaluate_shape_bilinear( xi, eta, phi, dphidxi, dphideta )
  real,               intent(in)  :: xi  !< The x position to evaluate
  real,               intent(in)  :: eta !< The z position to evaluate
  real, dimension(4), intent(out) :: phi !< The weights of the four corners at this point
  real, dimension(4), intent(out) :: dphidxi  !< The x-gradient of the weights of the four
                                         !! corners at this point
  real, dimension(4), intent(out) :: dphideta !< The z-gradient of the weights of the four
                                         !! corners at this point
 
  ! The shape functions within the parent element are defined as shown here:
  !
  !    (-1,1) 2 o------------o 1 (1,1)
  !             |            |
  !             |            |
  !             |            |
  !             |            |
  !   (-1,-1) 3 o------------o 4 (1,-1)
  !
 
  phi(1) = 0.25 * ( 1 + xi ) * ( 1 + eta )
  phi(2) = 0.25 * ( 1 - xi ) * ( 1 + eta )
  phi(3) = 0.25 * ( 1 - xi ) * ( 1 - eta )
  phi(4) = 0.25 * ( 1 + xi ) * ( 1 - eta )
 
  dphidxi(1) = 0.25 * ( 1 + eta )
  dphidxi(2) = - 0.25 * ( 1 + eta )
  dphidxi(3) = - 0.25 * ( 1 - eta )
  dphidxi(4) = 0.25 * ( 1 - eta )
 
  dphideta(1) = 0.25 * ( 1 + xi )
  dphideta(2) = 0.25 * ( 1 - xi )
  dphideta(3) = - 0.25 * ( 1 - xi )
  dphideta(4) = - 0.25 * ( 1 + xi )
 
end subroutine evaluate_shape_bilinear
 
 
! =============================================================================
!> Evaluation of the nine quadratic shape fn weights and their gradients at (xi,eta)
subroutine evaluate_shape_quadratic ( xi, eta, phi, dphidxi, dphideta )
 
  ! Arguments
  real,               intent(in)  :: xi  !< The x position to evaluate
  real,               intent(in)  :: eta !< The z position to evaluate
  real, dimension(9), intent(out) :: phi !< The weights of the 9 bilinear quadrature points
                                         !! at this point
  real, dimension(9), intent(out) :: dphidxi  !< The x-gradient of the weights of the 9 bilinear
                                         !! quadrature points corners at this point
  real, dimension(9), intent(out) :: dphideta !< The z-gradient of the weights of the 9 bilinear
                                         !! quadrature points corners at this point
 
  ! The quadratic shape functions within the parent element are defined as shown here:
  !
  !                 5 (0,1)
  !    (-1,1) 2 o------o------o 1 (1,1)
  !             |             |
  !             |   9 (0,0)   |
  !    (-1,0) 6 o      o      o 8 (1,0)
  !             |             |
  !             |             |
  !   (-1,-1) 3 o------o------o 4 (1,-1)
  !                 7 (0,-1)
  !
 
  phi(:)   = 0.0
  dphidxi(:)  = 0.0
  dphideta(:) = 0.0
 
  phi(1) = 0.25 * xi * ( 1 + xi ) * eta * ( 1 + eta )
  phi(2) = - 0.25 * xi * ( 1 - xi ) * eta * ( 1 + eta )
  phi(3) = 0.25 * xi * ( 1 - xi ) * eta * ( 1 - eta )
  phi(4) = - 0.25 * xi * ( 1 + xi ) * eta * ( 1 - eta )
  phi(5) = 0.5 * ( 1 + xi ) * ( 1 - xi ) * eta * ( 1 + eta )
  phi(6) = - 0.5 * xi * ( 1 - xi ) * ( 1 - eta ) * ( 1 + eta )
  phi(7) = - 0.5 * ( 1 - xi ) * ( 1 + xi ) * eta * ( 1 - eta )
  phi(8) = 0.5 * xi * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
  phi(9) = ( 1 - xi ) * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
 
  !dphidxi(1) = 0.25 * ( 1 + 2*xi ) * eta * ( 1 + eta )
  !dphidxi(2) = - 0.25 * ( 1 - 2*xi ) * eta * ( 1 + eta )
  !dphidxi(3) = 0.25 * ( 1 - 2*xi ) * eta * ( 1 - eta )
  !dphidxi(4) = - 0.25 * ( 1 + 2*xi ) * eta * ( 1 - eta )
  !dphidxi(5) = - xi * eta * ( 1 + eta )
  !dphidxi(6) = - 0.5 * ( 1 - 2*xi ) * ( 1 - eta ) * ( 1 + eta )
  !dphidxi(7) = xi * eta * ( 1 - eta )
  !dphidxi(8) = 0.5 * ( 1 + 2*xi ) * ( 1 - eta ) * ( 1 + eta )
  !dphidxi(9) = - 2 * xi * ( 1 - eta ) * ( 1 + eta )
 
  !dphideta(1) = 0.25 * xi * ( 1 + xi ) * ( 1 + 2*eta )
  !dphideta(2) = - 0.25 * xi * ( 1 - xi ) * ( 1 + 2*eta )
  !dphideta(3) = 0.25 * xi * ( 1 - xi ) * ( 1 - 2*eta )
  !dphideta(4) = - 0.25 * xi * ( 1 + xi ) * ( 1 - 2*eta )
  !dphideta(5) = 0.5 * ( 1 + xi ) * ( 1 - xi ) * ( 1 + 2*eta )
  !dphideta(6) = xi * ( 1 - xi ) * eta
  !dphideta(7) = - 0.5 * ( 1 - xi ) * ( 1 + xi ) * ( 1 - 2*eta )
  !dphideta(8) = - xi * ( 1 + xi ) * eta
  !dphideta(9) = - 2 * ( 1 - xi ) * ( 1 + xi ) * eta
 
end subroutine evaluate_shape_quadratic
! ==============================================================================
 
!>   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 Bode's rule quadrature to do the integrals.
subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, &
                                   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(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: t  !< Potential temperature of the layer [degC].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: s  !< Salinity of the layer [ppt].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_t !< Pressure atop the layer [Pa].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_b !< Pressure below the layer [Pa].
  real,                 intent(in)  :: alpha_ref !< A mean specific volume that is
                            !! subtracted out to reduce the magnitude of each of the
                            !! integrals [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
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(out) :: dza !< The change in the geopotential anomaly
                            !! across the layer [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(out) :: intp_dza !< The integral in pressure through the
                            !! layer of the geopotential anomaly relative to the anomaly
                            !! at the bottom of the layer [Pa m2 s-2].
  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
              optional, intent(out) :: 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 [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
              optional, intent(out) :: 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 [m2 s-2].
  integer,    optional, intent(in)  :: halo_size !< The width of halo points on which to calculate dza.
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(in)  :: bathyp !< The pressure at the bathymetry [Pa]
  real,       optional, intent(in)  :: dp_neglect !< A miniscule pressure change with
                                             !! the same units as p_t (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
! Bode'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), s5(5), p5(5), a5(5)
  real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1].
  real :: dp         ! The pressure change through a layer [Pa].
!  real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa].
  real :: hwght      ! A pressure-thickness below topography [Pa].
  real :: hl, hr     ! Pressure-thicknesses of the columns to the left and right [Pa].
  real :: idenom     ! The inverse of the denominator in the weights [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 [m2 s-2].
  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
 
  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) = p_b(i,j) - 0.25*real(n-1)*dp
    enddo
    call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode'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 Bode'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 wekghted.
      p5(1) = 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) - 0.25*dp
      enddo
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode'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 Bode'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 wekghted.
      p5(1) = 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) = p5(n-1) - 0.25*dp
      enddo
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode'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 Bode'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
 
!>   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 Bode'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, dza, &
                             intp_dza, intx_dza, inty_dza, useMassWghtInterp)
  type(hor_index_type), intent(in)  :: hi !< A horizontal index type structure.
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: t_t  !< Potential temperature at the top of the layer [degC].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: t_b  !< Potential temperature at the bottom of the layer [degC].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: s_t  !< Salinity at the top the layer [ppt].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: s_b  !< Salinity at the bottom the layer [ppt].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_t !< Pressure atop the layer [Pa].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_b !< Pressure below the layer [Pa].
  real,                 intent(in)  :: alpha_ref !< A mean specific volume that is
                            !! subtracted out to reduce the magnitude of each of the
                            !! integrals [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 (Pa?)
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: bathyp !< The pressure at the bathymetry [Pa]
  type(eos_type),       pointer     :: eos !< Equation of state structure
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(out) :: dza !< The change in the geopotential anomaly
                            !! across the layer [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(out) :: intp_dza !< The integral in pressure through the
                            !! layer of the geopotential anomaly relative to the anomaly
                            !! at the bottom of the layer [Pa m2 s-2].
  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
              optional, intent(out) :: 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 [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
              optional, intent(out) :: 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 [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
! Bode'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, dimension(5) :: t5, s5, p5, a5
  real, dimension(15) :: t15, s15, p15, a15
  real :: wt_t(5), wt_b(5)
  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 [Pa].
  real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa].
  real :: hwght      ! A pressure-thickness below topography [Pa].
  real :: hl, hr     ! Pressure-thicknesses of the columns to the left and right [Pa].
  real :: idenom     ! The inverse of the denominator in the weights [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 [m2 s-2].
  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
 
  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) = 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
    call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode'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 Bode'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 wekghted.
      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) = 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
 
    call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref)
 
    intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
    do m=2,4
      ! Use Bode'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 Bode'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 wekghted.
      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) = 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
 
    call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref)
 
    intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
    do m=2,4
      ! Use Bode'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 Bode'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
 
!> Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10
subroutine convert_temp_salt_for_teos10(T, S, press, G, kd, mask_z, EOS)
  use mom_grid, only : ocean_grid_type
 
  type(ocean_grid_type), intent(in)    :: g   !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                         intent(inout) :: t   !< Potential temperature referenced to the surface [degC]
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                         intent(inout) :: s   !< Salinity [ppt]
  real, dimension(:),    intent(in)    :: press !< Pressure at the top of the layer [Pa].
  type(eos_type),        pointer       :: eos !< Equation of state structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                         intent(in)    :: mask_z !< 3d mask regulating which points to convert.
  integer,               intent(in)    :: kd  !< The number of layers to work on
 
  integer :: i,j,k
  real :: gsw_sr_from_sp, gsw_ct_from_pt, gsw_sa_from_sp
  real :: p
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "convert_temp_salt_to_TEOS10 called with an unassociated EOS_type EOS.")
 
  if ((eos%form_of_EOS /= eos_teos10) .and. (eos%form_of_EOS /= eos_nemo)) return
 
  do k=1,kd ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if (mask_z(i,j,k) >= 1.0) then
     s(i,j,k) = gsw_sr_from_sp(s(i,j,k))
!     p=press(k)/10000. !convert pascal to dbar
!     S(i,j,k) = gsw_sa_from_sp(S(i,j,k),p,G%geoLonT(i,j),G%geoLatT(i,j))
     t(i,j,k) = gsw_ct_from_pt(s(i,j,k),t(i,j,k))
    endif
  enddo ; enddo ; enddo
end subroutine convert_temp_salt_for_teos10
 
!> Extractor routine for the EOS type if the members need to be accessed outside this module
subroutine extract_member_eos(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, &
                              Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
  type(eos_type),    pointer     :: eos !< Equation of state structure
  integer, optional, intent(out) :: form_of_eos !< A coded integer indicating the equation of state to use.
  integer, optional, intent(out) :: form_of_tfreeze !< A coded integer indicating the expression for
                                       !! the potential temperature of the freezing point.
  logical, optional, intent(out) :: eos_quadrature !< If true, always use the generic (quadrature)
                                       !! code for the integrals of density.
  logical, optional, intent(out) :: compressible !< If true, in situ density is a function of pressure.
  real   , optional, intent(out) :: rho_t0_s0 !< Density at T=0 degC and S=0 ppt [kg m-3]
  real   , optional, intent(out) :: drho_dt   !< Partial derivative of density with temperature
                                              !! in [kg m-3 degC-1]
  real   , optional, intent(out) :: drho_ds   !< Partial derivative of density with salinity
                                              !! in [kg m-3 ppt-1]
  real   , optional, intent(out) :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC].
  real   , optional, intent(out) :: dtfr_ds   !< The derivative of freezing point with salinity
                                              !! [degC PSU-1].
  real   , optional, intent(out) :: dtfr_dp   !< The derivative of freezing point with pressure
                                              !! [degC Pa-1].
 
  if (present(form_of_eos    ))  form_of_eos     = eos%form_of_EOS
  if (present(form_of_tfreeze))  form_of_tfreeze = eos%form_of_TFreeze
  if (present(eos_quadrature ))  eos_quadrature  = eos%EOS_quadrature
  if (present(compressible   ))  compressible    = eos%Compressible
  if (present(rho_t0_s0      ))  rho_t0_s0       = eos%Rho_T0_S0
  if (present(drho_dt        ))  drho_dt         = eos%drho_dT
  if (present(drho_ds        ))  drho_ds         = eos%dRho_dS
  if (present(tfr_s0_p0      ))  tfr_s0_p0       = eos%TFr_S0_P0
  if (present(dtfr_ds        ))  dtfr_ds         = eos%dTFr_dS
  if (present(dtfr_dp        ))  dtfr_dp         = eos%dTFr_dp
 
end subroutine extract_member_eos
 
end module mom_eos
 
!> \namespace mom_eos
!!
!! The MOM_EOS module is a wrapper for various equations of state (e.g. Linear,
!! Wright, UNESCO) and provides a uniform interface to the rest of the model
!! independent of which equation of state is being used.