MOM_PressureForce_analytic_FV.F90

!> Analytically integrated finite volume pressure gradient
module mom_pressureforce_afv
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_diag_mediator, only : post_data, register_diag_field
use mom_diag_mediator, only : safe_alloc_ptr, diag_ctrl, time_type
use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_grid, only : ocean_grid_type
use mom_pressureforce_mont, only : set_pbce_bouss, set_pbce_nonbouss
use mom_tidal_forcing, only : calc_tidal_forcing, tidal_forcing_cs
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs
use mom_eos, only : int_density_dz, int_specific_vol_dp
use mom_eos, only : int_density_dz_generic_plm, int_density_dz_generic_ppm
use mom_eos, only : int_spec_vol_dp_generic_plm
use mom_eos, only : int_density_dz_generic, int_spec_vol_dp_generic
use mom_ale, only : pressure_gradient_plm, pressure_gradient_ppm, ale_cs
 
implicit none ; private
 
#include <MOM_memory.h>
 
public pressureforce_afv, pressureforce_afv_init, pressureforce_afv_end
public pressureforce_afv_bouss, pressureforce_afv_nonbouss
 
! 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.
 
!> Finite volume pressure gradient control structure
type, public :: pressureforce_afv_cs ; private
  logical :: tides          !< If true, apply tidal momentum forcing.
  real    :: rho0           !< The density used in the Boussinesq
                            !! approximation [kg m-3].
  real    :: gfs_scale      !< A scaling of the surface pressure gradients to
                            !! allow the use of a reduced gravity model [nondim].
  type(time_type), pointer :: time !< A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                            !! timing of diagnostic output.
  logical :: usemasswghtinterp !< Use mass weighting in T/S interpolation
  logical :: boundary_extrap !< Indicate whether high-order boundary
                            !! extrapolation should be used within boundary cells
 
  logical :: reconstruct    !< If true, polynomial profiles of T & S will be
                            !! reconstructed and used in the integrals for the
                            !! finite volume pressure gradient calculation.
                            !! The default depends on whether regridding is being used.
 
  integer :: recon_scheme   !< Order of the polynomial of the reconstruction of T & S
                            !! for the finite volume pressure gradient calculation.
                            !! By the default (1) is for a piecewise linear method
 
  integer :: id_e_tidal = -1 !< Diagnostic identifier
  type(tidal_forcing_cs), pointer :: tides_csp => null() !< Tides control structure
end type pressureforce_afv_cs
 
contains
 
!> Thin interface between the model and the Boussinesq and non-Boussinesq
!! pressure force routines.
subroutine pressureforce_afv(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
  type(ocean_grid_type),                     intent(in)    :: g   !< Ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv  !< Vertical grid structure
  type(unit_scale_type),                     intent(in)    :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h   !< Layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                     intent(inout) :: tv  !< Thermodynamic variables
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out)   :: pfu !< Zonal acceleration [m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: pfv !< Meridional acceleration [m s-2]
  type(pressureforce_afv_cs),                pointer       :: cs  !< Finite volume PGF control structure
  type(ale_cs),                              pointer       :: ale_csp !< ALE control structure
  real, dimension(:,:),                      optional, pointer :: p_atm !< The pressure at the ice-ocean
                                                           !! or atmosphere-ocean interface [Pa].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  optional, intent(out) :: pbce !< The baroclinic pressure
                                                           !! anomaly in each layer due to eta anomalies
                                                           !! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
  real, dimension(SZI_(G),SZJ_(G)),          optional, intent(out) :: eta !< The bottom mass used to
                                                           !! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
                                                           !! contributions or compressibility compensation.
 
  if (gv%Boussinesq) then
    call pressureforce_afv_bouss(h, tv, pfu, pfv, g, gv, us, cs, ale_csp, p_atm, pbce, eta)
  else
    call pressureforce_afv_nonbouss(h, tv, pfu, pfv, g, gv, us, cs, ale_csp, p_atm, pbce, eta)
  endif
 
end subroutine pressureforce_afv
 
!> \brief Non-Boussinesq analytically-integrated finite volume form of pressure gradient
!!
!! Determines the acceleration due to hydrostatic pressure forces, using
!! the analytic finite volume form of the Pressure gradient, and does not
!! make the Boussinesq approximation.
!!
!! To work, the following fields must be set outside of the usual (is:ie,js:je)
!! range before this subroutine is called:
!!   h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
subroutine pressureforce_afv_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
  type(ocean_grid_type),                     intent(in)  :: g   !< Ocean grid structure
  type(verticalgrid_type),                   intent(in)  :: gv  !< Vertical grid structure
  type(unit_scale_type),                     intent(in)  :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)  :: h   !< Layer thickness [H ~> kg/m2]
  type(thermo_var_ptrs),                     intent(in)  :: tv  !< Thermodynamic variables
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: pfu !< Zonal acceleration [m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration [m s-2]
  type(pressureforce_afv_cs),                pointer     :: cs  !< Finite volume PGF control structure
  type(ale_cs),                              pointer     :: ale_csp !< ALE control structure
  real, dimension(:,:),                      optional, pointer :: p_atm !< The pressure at the ice-ocean
                                                           !! or atmosphere-ocean interface [Pa].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  optional, intent(out) :: pbce !< The baroclinic pressure
                                                           !! anomaly in each layer due to eta anomalies
                                                           !! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
  real, dimension(SZI_(G),SZJ_(G)),          optional, intent(out) :: eta !< The bottom mass used to
                                                           !! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
                                                           !! contributions or compressibility compensation.
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure [Pa].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
    t_tmp, &    ! Temporary array of temperatures where layers that are lighter
                ! than the mixed layer have the mixed layer's properties [degC].
    s_tmp       ! Temporary array of salinities where layers that are lighter
                ! than the mixed layer have the mixed layer's properties [ppt].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    s_t, &      ! Top and bottom edge values for linear reconstructions
    s_b, &      ! of salinity within each layer [ppt].
    t_t, &      ! Top and bottom edge values for linear reconstructions
    t_b         ! of temperature within each layer [degC].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: &
    dza, &      ! The change in geopotential anomaly between the top and bottom
                ! of a layer [m2 s-2].
    intp_dza    ! The vertical integral in depth of the pressure anomaly less
                ! the pressure anomaly at the top of the layer [Pa m2 s-2].
  real, dimension(SZI_(G),SZJ_(G))  :: &
    dp, &       ! The (positive) change in pressure across a layer [Pa].
    ssh, &      ! The sea surface height anomaly, in depth units [Z ~> m].
    e_tidal, &  ! The bottom geopotential anomaly due to tidal forces from
                ! astronomical sources and self-attraction and loading [Z ~> m].
    dm, &       ! The barotropic adjustment to the Montgomery potential to
                ! account for a reduced gravity model [m2 s-2].
    za          ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the
                ! interface atop a layer [m2 s-2].
 
  real, dimension(SZI_(G)) :: rho_cv_bl !  The coordinate potential density in the deepest variable
                ! density near-surface layer [kg m-3].
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    intx_za     ! The zonal integral of the geopotential anomaly along the
                ! interface below a layer, divided by the grid spacing [m2 s-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
    intx_dza    ! The change in intx_za through a layer [m2 s-2].
  real, dimension(SZI_(G),SZJB_(G)) :: &
    inty_za     ! The meridional integral of the geopotential anomaly along the
                ! interface below a layer, divided by the grid spacing [m2 s-2].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
    inty_dza    ! The change in inty_za through a layer [m2 s-2].
  real :: p_ref(szi_(g))     !   The pressure used to calculate the coordinate
                             ! density, [Pa] (usually 2e7 Pa = 2000 dbar).
 
  real :: dp_neglect         ! A thickness that is so small it is usually lost
                             ! in roundoff and can be neglected [Pa].
  real :: g_earth_z          ! A scaled version of g_Earth [m2 Z-1 s-2 ~> m s-2].
  real :: i_gearth           ! The inverse of g_Earth_z [s2 Z m-2 ~> s2 m-1]
  real :: alpha_anom         ! The in-situ specific volume, averaged over a
                             ! layer, less alpha_ref [m3 kg-1].
  logical :: use_p_atm       ! If true, use the atmospheric pressure.
  logical :: use_ale         ! If true, use an ALE pressure reconstruction.
  logical :: use_eos    ! If true, density is calculated from T & S using an
                        ! equation of state.
  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
 
  real :: alpha_ref ! A reference specific volume [m3 kg-1], that is used
                    ! to reduce the impact of truncation errors.
  real :: rho_in_situ(szi_(g)) ! The in situ density [kg m-3].
  real :: pa_to_h   ! A factor to convert from Pa to the thicknesss units (H).
!  real :: oneatm = 101325.0  ! 1 atm in [Pa] = [kg m-1 s-2]
  real, parameter :: c1_6 = 1.0/6.0
  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
  integer :: i, j, k
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  nkmb=gv%nk_rho_varies
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
 
  if (.not.associated(cs)) call mom_error(fatal, &
       "MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.")
 
  use_p_atm = .false.
  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
  use_eos = associated(tv%eqn_of_state)
  use_ale = .false.
  if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
 
  dp_neglect = gv%H_to_Pa * gv%H_subroundoff
  alpha_ref = 1.0/cs%Rho0
  g_earth_z = us%L_T_to_m_s**2 * gv%g_Earth
  i_gearth = 1.0 / g_earth_z
 
  if (use_p_atm) then
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      p(i,j,1) = p_atm(i,j)
    enddo ; enddo
  else
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      p(i,j,1) = 0.0 ! or oneatm
    enddo ; enddo
  endif
  !$OMP parallel do default(shared)
  do j=jsq,jeq+1 ; do k=2,nz+1 ; do i=isq,ieq+1
    p(i,j,k) = p(i,j,k-1) + gv%H_to_Pa * h(i,j,k-1)
  enddo ; enddo ; enddo
 
  if (use_eos) then
  !   With a bulk mixed layer, replace the T & S of any layers that are
  ! lighter than the the buffer layer with the properties of the buffer
  ! layer.  These layers will be massless anyway, and it avoids any
  ! formal calculations with hydrostatically unstable profiles.
    if (nkmb>0) then
      tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
      tv_tmp%eqn_of_state => tv%eqn_of_state
      do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
      !$OMP parallel do default(shared) private(Rho_cv_BL)
      do j=jsq,jeq+1
        do k=1,nkmb ; do i=isq,ieq+1
          tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
        enddo ; enddo
        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
                        rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state)
        do k=nkmb+1,nz ; do i=isq,ieq+1
          if (gv%Rlay(k) < rho_cv_bl(i)) then
            tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
          else
            tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
          endif
        enddo ; enddo
      enddo
    else
      tv_tmp%T => tv%T ; tv_tmp%S => tv%S
      tv_tmp%eqn_of_state => tv%eqn_of_state
    endif
  endif
 
  ! If regridding is activated, do a linear reconstruction of salinity
  ! and temperature across each layer. The subscripts 't' and 'b' refer
  ! to top and bottom values within each layer (these are the only degrees
  ! of freedeom needed to know the linear profile).
  if ( use_ale ) then
    if ( cs%Recon_Scheme == 1 ) then
      call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
    elseif ( cs%Recon_Scheme == 2) then
      call pressure_gradient_ppm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
    endif
  endif
 
  !$OMP parallel do default(shared) private(alpha_anom,dp)
  do k=1,nz
    ! Calculate 4 integrals through the layer that are required in the
    ! subsequent calculation.
    if (use_eos) then
      if ( use_ale ) then
        if ( cs%Recon_Scheme == 1 ) then
          call int_spec_vol_dp_generic_plm( t_t(:,:,k), t_b(:,:,k), &
                    s_t(:,:,k), s_b(:,:,k), p(:,:,k), p(:,:,k+1), &
                    alpha_ref, dp_neglect, p(:,:,nz+1), g%HI, &
                    tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
                    intx_dza(:,:,k), inty_dza(:,:,k), &
                    usemasswghtinterp = cs%useMassWghtInterp)
          i=k
        elseif ( cs%Recon_Scheme == 2 ) then
          call mom_error(fatal, "PressureForce_AFV_nonBouss: "//&
                         "int_spec_vol_dp_generic_ppm does not exist yet.")
        !  call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
        !            tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), &
        !            alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
        !            intx_dza(:,:,k), inty_dza(:,:,k))
        endif
      else
        call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,k), &
                               p(:,:,k+1), alpha_ref, g%HI, tv%eqn_of_state, &
                               dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
                               inty_dza(:,:,k), bathyp=p(:,:,nz+1), dp_tiny=dp_neglect, &
                               usemasswghtinterp = cs%useMassWghtInterp)
      endif
    else
      alpha_anom = 1.0/gv%Rlay(k) - alpha_ref
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        dp(i,j) = gv%H_to_Pa * h(i,j,k)
        dza(i,j,k) = alpha_anom * dp(i,j)
        intp_dza(i,j,k) = 0.5 * alpha_anom * dp(i,j)**2
      enddo ; enddo
      do j=js,je ; do i=isq,ieq
        intx_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i+1,j))
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        inty_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i,j+1))
      enddo ; enddo
    endif
  enddo
 
  !   The bottom geopotential anomaly is calculated first so that the increments
  ! to the geopotential anomalies can be reused.  Alternately, the surface
  ! geopotential could be calculated directly with separate calls to
  ! int_specific_vol_dp with alpha_ref=0, and the anomalies used going
  ! downward, which would relieve the need for dza, intp_dza, intx_dza, and
  ! inty_dza to be 3-D arrays.
 
  ! Sum vertically to determine the surface geopotential anomaly.
  !$OMP parallel do default(shared)
  do j=jsq,jeq+1
    do i=isq,ieq+1
      za(i,j) = alpha_ref*p(i,j,nz+1) - g_earth_z*g%bathyT(i,j)
    enddo
    do k=nz,1,-1 ; do i=isq,ieq+1
      za(i,j) = za(i,j) + dza(i,j,k)
    enddo ; enddo
  enddo
 
  if (cs%tides) then
    ! Find and add the tidal geopotential anomaly.
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      ssh(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * i_gearth
    enddo ; enddo
    call calc_tidal_forcing(cs%Time, ssh, e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      za(i,j) = za(i,j) - g_earth_z * e_tidal(i,j)
    enddo ; enddo
  endif
 
  if (cs%GFS_scale < 1.0) then
    ! Adjust the Montgomery potential to make this a reduced gravity model.
    if (use_eos) then
      !$OMP parallel do default(shared) private(rho_in_situ)
      do j=jsq,jeq+1
        call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), &
                               rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
 
        do i=isq,ieq+1
          dm(i,j) = (cs%GFS_scale - 1.0) * &
            (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
        enddo
      enddo
    else
      !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        dm(i,j) = (cs%GFS_scale - 1.0) * &
          (p(i,j,1)*(1.0/gv%Rlay(1) - alpha_ref) + za(i,j))
      enddo ; enddo
    endif
!  else
!    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; dM(i,j) = 0.0 ; enddo ; enddo
  endif
 
  !   This order of integrating upward and then downward again is necessary with
  ! a nonlinear equation of state, so that the surface geopotentials will go
  ! linearly between the values at thickness points, but the bottom
  ! geopotentials will not now be linear at the sub-grid-scale.  Doing this
  ! ensures no motion with flat isopycnals, even with a nonlinear equation of state.
  !$OMP parallel do default(shared)
  do j=js,je ; do i=isq,ieq
    intx_za(i,j) = 0.5*(za(i,j) + za(i+1,j))
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=jsq,jeq ; do i=is,ie
    inty_za(i,j) = 0.5*(za(i,j) + za(i,j+1))
  enddo ; enddo
  do k=1,nz
    ! These expressions for the acceleration have been carefully checked in
    ! a set of idealized cases, and should be bug-free.
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      dp(i,j) = gv%H_to_Pa*h(i,j,k)
      za(i,j) = za(i,j) - dza(i,j,k)
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=js,je ; do i=isq,ieq
      intx_za(i,j) = intx_za(i,j) - intx_dza(i,j,k)
      pfu(i,j,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
                   (za(i+1,j)*dp(i+1,j) + intp_dza(i+1,j,k))) + &
                   ((dp(i+1,j) - dp(i,j)) * intx_za(i,j) - &
                   (p(i+1,j,k) - p(i,j,k)) * intx_dza(i,j,k))) * &
                   (2.0*g%IdxCu(i,j) / ((dp(i,j) + dp(i+1,j)) + &
                   dp_neglect))
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=jsq,jeq ; do i=is,ie
      inty_za(i,j) = inty_za(i,j) - inty_dza(i,j,k)
      pfv(i,j,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
                   (za(i,j+1)*dp(i,j+1) + intp_dza(i,j+1,k))) + &
                   ((dp(i,j+1) - dp(i,j)) * inty_za(i,j) - &
                   (p(i,j+1,k) - p(i,j,k)) * inty_dza(i,j,k))) * &
                   (2.0*g%IdyCv(i,j) / ((dp(i,j) + dp(i,j+1)) + &
                   dp_neglect))
    enddo ; enddo
 
    if (cs%GFS_scale < 1.0) then
      ! Adjust the Montgomery potential to make this a reduced gravity model.
      !$OMP parallel do default(shared)
      do j=js,je ; do i=isq,ieq
        pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
      enddo ; enddo
      !$OMP parallel do default(shared)
      do j=jsq,jeq ; do i=is,ie
        pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
      enddo ; enddo
    endif
  enddo
 
  if (present(pbce)) then
    call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce)
  endif
 
  if (present(eta)) then
    pa_to_h = 1.0 / gv%H_to_Pa
    if (use_p_atm) then
      !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h ! eta has the same units as h.
      enddo ; enddo
    else
      !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        eta(i,j) = p(i,j,nz+1)*pa_to_h ! eta has the same units as h.
      enddo ; enddo
    endif
  endif
 
  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
 
end subroutine pressureforce_afv_nonbouss
 
!> \brief Boussinesq analytically-integrated finite volume form of pressure gradient
!!
!! Determines the acceleration due to hydrostatic pressure forces, using
!! the finite volume form of the terms and analytic integrals in depth.
!!
!! To work, the following fields must be set outside of the usual (is:ie,js:je)
!! range before this subroutine is called:
!!   h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
subroutine pressureforce_afv_bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
  type(ocean_grid_type),                     intent(in)  :: g   !< Ocean grid structure
  type(verticalgrid_type),                   intent(in)  :: gv  !< Vertical grid structure
  type(unit_scale_type),                     intent(in)  :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)  :: h   !< Layer thickness [H ~> m]
  type(thermo_var_ptrs),                     intent(in)  :: tv  !< Thermodynamic variables
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: pfu !< Zonal acceleration [m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration [m s-2]
  type(pressureforce_afv_cs),                pointer     :: cs  !< Finite volume PGF control structure
  type(ale_cs),                              pointer     :: ale_csp !< ALE control structure
  real, dimension(:,:),                      optional, pointer :: p_atm !< The pressure at the ice-ocean
                                                         !! or atmosphere-ocean interface [Pa].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  optional, intent(out) :: pbce !< The baroclinic pressure
                                                         !! anomaly in each layer due to eta anomalies
                                                         !! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
  real, dimension(SZI_(G),SZJ_(G)),          optional, intent(out) :: eta !< The bottom mass used to
                                                         !! calculate PFu and PFv [H ~> m or kg m-2], with any
                                                         !! tidal contributions or compressibility compensation.
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in depth units [Z ~> m].
  real, dimension(SZI_(G),SZJ_(G))  :: &
    e_tidal, &  ! The bottom geopotential anomaly due to tidal forces from
                ! astronomical sources and self-attraction and loading [Z ~> m].
    dm          ! The barotropic adjustment to the Montgomery potential to
                ! account for a reduced gravity model [m2 s-2].
  real, dimension(SZI_(G)) :: &
    rho_cv_bl   !   The coordinate potential density in the deepest variable
                ! density near-surface layer [kg m-3].
  real, dimension(SZI_(G),SZJ_(G)) :: &
    dz, &       ! The change in geopotential thickness through a layer [m2 s-2].
    pa, &       ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the
                ! the interface atop a layer [Pa].
    dpa, &      ! The change in pressure anomaly between the top and bottom
                ! of a layer [Pa].
    intz_dpa    ! The vertical integral in depth of the pressure anomaly less the
                ! pressure anomaly at the top of the layer [H Pa ~> m Pa or kg m-2 Pa].
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    intx_pa, &  ! The zonal integral of the pressure anomaly along the interface
                ! atop a layer, divided by the grid spacing [Pa].
    intx_dpa    ! The change in intx_pa through a layer [Pa].
  real, dimension(SZI_(G),SZJB_(G)) :: &
    inty_pa, &  ! The meridional integral of the pressure anomaly along the
                ! interface atop a layer, divided by the grid spacing [Pa].
    inty_dpa    ! The change in inty_pa through a layer [Pa].
 
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
    t_tmp, &    ! Temporary array of temperatures where layers that are lighter
                ! than the mixed layer have the mixed layer's properties [degC].
    s_tmp       ! Temporary array of salinities where layers that are lighter
                ! than the mixed layer have the mixed layer's properties [ppt].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    s_t, s_b, t_t, t_b ! Top and bottom edge values for linear reconstructions
                       ! of salinity and temperature within each layer.
  real :: rho_in_situ(szi_(g)) ! The in situ density [kg m-3].
  real :: p_ref(szi_(g))     !   The pressure used to calculate the coordinate
                             ! density, [Pa] (usually 2e7 Pa = 2000 dbar).
  real :: p0(szi_(g))        ! An array of zeros to use for pressure [Pa].
  real :: h_neglect          ! A thickness that is so small it is usually lost
                             ! in roundoff and can be neglected [H ~> m].
  real :: g_earth_z          ! A scaled version of g_Earth [m2 Z-1 s-2 ~> m s-2].
  real :: i_rho0             ! 1/Rho0 [m3 kg-1].
  real :: g_rho0             ! G_Earth / Rho0 in [m5 Z-1 s-2 kg-1 ~> m4 s-2 kg-1].
  real :: rho_ref            ! The reference density [kg m-3].
  real :: dz_neglect         ! A minimal thickness [Z ~> m], like e.
  logical :: use_p_atm       ! If true, use the atmospheric pressure.
  logical :: use_ale         ! If true, use an ALE pressure reconstruction.
  logical :: use_eos         ! If true, density is calculated from T & S using an equation of state.
  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
 
  real, parameter :: c1_6 = 1.0/6.0
  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
  integer :: i, j, k
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  nkmb=gv%nk_rho_varies
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
 
  if (.not.associated(cs)) call mom_error(fatal, &
       "MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.")
 
  use_p_atm = .false.
  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
  use_eos = associated(tv%eqn_of_state)
  do i=isq,ieq+1 ; p0(i) = 0.0 ; enddo
  use_ale = .false.
  if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
 
  h_neglect = gv%H_subroundoff
  dz_neglect = gv%H_subroundoff * gv%H_to_Z
  i_rho0 = 1.0/gv%Rho0
  g_earth_z = us%L_T_to_m_s**2 * gv%g_Earth
  g_rho0 = g_earth_z/gv%Rho0
  rho_ref = cs%Rho0
 
  if (cs%tides) then
    !   Determine the surface height anomaly for calculating self attraction
    ! and loading.  This should really be based on bottom pressure anomalies,
    ! but that is not yet implemented, and the current form is correct for
    ! barotropic tides.
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1
      do i=isq,ieq+1
        e(i,j,1) = -g%bathyT(i,j)
      enddo
      do k=1,nz ; do i=isq,ieq+1
        e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
      enddo ; enddo
    enddo
    call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
  endif
 
!    Here layer interface heights, e, are calculated.
  if (cs%tides) then
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
    enddo ; enddo
  else
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      e(i,j,nz+1) = -g%bathyT(i,j)
    enddo ; enddo
  endif
  !$OMP parallel do default(shared)
  do j=jsq,jeq+1; do k=nz,1,-1 ; do i=isq,ieq+1
    e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
  enddo ; enddo ; enddo
 
  if (use_eos) then
! With a bulk mixed layer, replace the T & S of any layers that are
! lighter than the the buffer layer with the properties of the buffer
! layer.  These layers will be massless anyway, and it avoids any
! formal calculations with hydrostatically unstable profiles.
 
    if (nkmb>0) then
      tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
      tv_tmp%eqn_of_state => tv%eqn_of_state
 
      do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
     !$OMP parallel do default(shared) private(Rho_cv_BL)
      do j=jsq,jeq+1
        do k=1,nkmb ; do i=isq,ieq+1
          tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
        enddo ; enddo
        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
                        rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state)
 
        do k=nkmb+1,nz ; do i=isq,ieq+1
          if (gv%Rlay(k) < rho_cv_bl(i)) then
            tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
          else
            tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
          endif
        enddo ; enddo
      enddo
    else
      tv_tmp%T => tv%T ; tv_tmp%S => tv%S
      tv_tmp%eqn_of_state => tv%eqn_of_state
    endif
  endif
 
  if (cs%GFS_scale < 1.0) then
    ! Adjust the Montgomery potential to make this a reduced gravity model.
    if (use_eos) then
      !$OMP parallel do default(shared)
      do j=jsq,jeq+1
        if (use_p_atm) then
          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), &
                                 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
        else
          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, &
                                 rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
        endif
        do i=isq,ieq+1
          dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * rho_in_situ(i)) * e(i,j,1)
        enddo
      enddo
    else
      !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * gv%Rlay(1)) * e(i,j,1)
      enddo ; enddo
    endif
  endif
  ! I have checked that rho_0 drops out and that the 1-layer case is right. RWH.
 
  ! If regridding is activated, do a linear reconstruction of salinity
  ! and temperature across each layer. The subscripts 't' and 'b' refer
  ! to top and bottom values within each layer (these are the only degrees
  ! of freedeom needed to know the linear profile).
  if ( use_ale ) then
    if ( cs%Recon_Scheme == 1 ) then
      call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
    elseif ( cs%Recon_Scheme == 2 ) then
      call pressure_gradient_ppm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
    endif
  endif
 
  ! Set the surface boundary conditions on pressure anomaly and its horizontal
  ! integrals, assuming that the surface pressure anomaly varies linearly
  ! in x and y.
  if (use_p_atm) then
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      pa(i,j) = (rho_ref*g_earth_z)*e(i,j,1) + p_atm(i,j)
    enddo ; enddo
  else
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      pa(i,j) = (rho_ref*g_earth_z)*e(i,j,1)
    enddo ; enddo
  endif
  !$OMP parallel do default(shared)
  do j=js,je ; do i=isq,ieq
    intx_pa(i,j) = 0.5*(pa(i,j) + pa(i+1,j))
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=jsq,jeq ; do i=is,ie
    inty_pa(i,j) = 0.5*(pa(i,j) + pa(i,j+1))
  enddo ; enddo
 
  do k=1,nz
    ! Calculate 4 integrals through the layer that are required in the
    ! subsequent calculation.
 
    if (use_eos) then
      ! The following routine computes the integrals that are needed to
      ! calculate the pressure gradient force. Linear profiles for T and S are
      ! assumed when regridding is activated. Otherwise, the previous version
      ! is used, whereby densities within each layer are constant no matter
      ! where the layers are located.
      if ( use_ale ) then
        if ( cs%Recon_Scheme == 1 ) then
          call int_density_dz_generic_plm( t_t(:,:,k), t_b(:,:,k), &
                    s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
                    rho_ref, cs%Rho0, g_earth_z, &
                    dz_neglect, g%bathyT, g%HI, g%HI, &
                    tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, &
                    usemasswghtinterp = cs%useMassWghtInterp)
        elseif ( cs%Recon_Scheme == 2 ) then
          call int_density_dz_generic_ppm( tv%T(:,:,k), t_t(:,:,k), t_b(:,:,k), &
                    tv%S(:,:,k), s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
                    rho_ref, cs%Rho0, g_earth_z, &
                    g%HI, g%HI, tv%eqn_of_state, dpa, intz_dpa, &
                    intx_dpa, inty_dpa)
        endif
      else
        call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,k), e(:,:,k+1), &
                  rho_ref, cs%Rho0, g_earth_z, g%HI, g%HI, tv%eqn_of_state, &
                  dpa, intz_dpa, intx_dpa, inty_dpa, &
                  g%bathyT, dz_neglect, cs%useMassWghtInterp)
      endif
      !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        intz_dpa(i,j) = intz_dpa(i,j)*gv%Z_to_H
      enddo ; enddo
    else
      !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        dz(i,j) = g_earth_z * gv%H_to_Z*h(i,j,k)
        dpa(i,j) = (gv%Rlay(k) - rho_ref)*dz(i,j)
        intz_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref)*dz(i,j)*h(i,j,k)
      enddo ; enddo
      !$OMP parallel do default(shared)
      do j=js,je ; do i=isq,ieq
        intx_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * (dz(i,j)+dz(i+1,j))
      enddo ; enddo
      !$OMP parallel do default(shared)
      do j=jsq,jeq ; do i=is,ie
        inty_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * (dz(i,j)+dz(i,j+1))
      enddo ; enddo
    endif
 
    ! Compute pressure gradient in x direction
    !$OMP parallel do default(shared)
    do j=js,je ; do i=isq,ieq
      pfu(i,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
                   (pa(i+1,j)*h(i+1,j,k) + intz_dpa(i+1,j))) + &
                   ((h(i+1,j,k) - h(i,j,k)) * intx_pa(i,j) - &
                   (e(i+1,j,k+1) - e(i,j,k+1)) * intx_dpa(i,j) * gv%Z_to_H)) * &
                   ((2.0*i_rho0*g%IdxCu(i,j)) / &
                   ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
      intx_pa(i,j) = intx_pa(i,j) + intx_dpa(i,j)
    enddo ; enddo
    ! Compute pressure gradient in y direction
    !$OMP parallel do default(shared)
    do j=jsq,jeq ; do i=is,ie
      pfv(i,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
                   (pa(i,j+1)*h(i,j+1,k) + intz_dpa(i,j+1))) + &
                   ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,j) - &
                   (e(i,j+1,k+1) - e(i,j,k+1)) * inty_dpa(i,j) * gv%Z_to_H)) * &
                   ((2.0*i_rho0*g%IdyCv(i,j)) / &
                   ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
      inty_pa(i,j) = inty_pa(i,j) + inty_dpa(i,j)
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      pa(i,j) = pa(i,j) + dpa(i,j)
    enddo ; enddo
  enddo
 
  if (cs%GFS_scale < 1.0) then
    do k=1,nz
      !$OMP parallel do default(shared)
      do j=js,je ; do i=isq,ieq
        pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
      enddo ; enddo
      !$OMP parallel do default(shared)
      do j=jsq,jeq ; do i=is,ie
        pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
      enddo ; enddo
    enddo
  endif
 
  if (present(pbce)) then
    call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce)
  endif
 
  if (present(eta)) then
    if (cs%tides) then
    ! eta is the sea surface height relative to a time-invariant geoid, for
    ! comparison with what is used for eta in btstep.  See how e was calculated
    ! about 200 lines above.
    !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
      enddo ; enddo
    else
    !$OMP parallel do default(shared)
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        eta(i,j) = e(i,j,1)*gv%Z_to_H
      enddo ; enddo
    endif
  endif
 
  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
 
end subroutine pressureforce_afv_bouss
 
!> Initializes the finite volume pressure gradient control structure
subroutine pressureforce_afv_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
  type(time_type), target,    intent(in)    :: time !< Current model time
  type(ocean_grid_type),      intent(in)    :: g  !< Ocean grid structure
  type(verticalgrid_type),    intent(in)    :: gv !< Vertical grid structure
  type(unit_scale_type),      intent(in)    :: us !< A dimensional unit scaling type
  type(param_file_type),      intent(in)    :: param_file !< Parameter file handles
  type(diag_ctrl), target,    intent(inout) :: diag !< Diagnostics control structure
  type(pressureforce_afv_cs), pointer       :: cs !< Finite volume PGF control structure
  type(tidal_forcing_cs), optional, pointer :: tides_csp !< Tides control structure
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl  ! This module's name.
  logical :: use_ale
 
  if (associated(cs)) then
    call mom_error(warning, "PressureForce_init called with an associated "// &
                            "control structure.")
    return
  else ; allocate(cs) ; endif
 
  cs%diag => diag ; cs%Time => time
  if (present(tides_csp)) then
    if (associated(tides_csp)) cs%tides_CSp => tides_csp
  endif
 
  mdl = "MOM_PressureForce_AFV"
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
                 "The mean ocean density used with BOUSSINESQ true to "//&
                 "calculate accelerations and the mass for conservation "//&
                 "properties, or with BOUSSINSEQ false to convert some "//&
                 "parameters from vertical units of m to kg m-2.", &
                 units="kg m-3", default=1035.0)
  call get_param(param_file, mdl, "TIDES", cs%tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)
  call get_param(param_file, "MOM", "USE_REGRIDDING", use_ale, &
                 "If True, use the ALE algorithm (regridding/remapping). "//&
                 "If False, use the layered isopycnal algorithm.", default=.false. )
  call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", cs%useMassWghtInterp, &
                 "If true, use mass weighting when interpolating T/S for "//&
                 "integrals near the bathymetry in AFV pressure gradient "//&
                 "calculations.", default=.false.)
  call get_param(param_file, mdl, "RECONSTRUCT_FOR_PRESSURE", cs%reconstruct, &
                 "If True, use vertical reconstruction of T & S within "//&
                 "the integrals of the FV pressure gradient calculation. "//&
                 "If False, use the constant-by-layer algorithm. "//&
                 "The default is set by USE_REGRIDDING.", &
                 default=use_ale )
  call get_param(param_file, mdl, "PRESSURE_RECONSTRUCTION_SCHEME", cs%Recon_Scheme, &
                 "Order of vertical reconstruction of T/S to use in the "//&
                 "integrals within the FV pressure gradient calculation.\n"//&
                 " 0: PCM or no reconstruction.\n"//&
                 " 1: PLM reconstruction.\n"//&
                 " 2: PPM reconstruction.", default=1)
  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION_PRESSURE", cs%boundary_extrap, &
                 "If true, the reconstruction of T & S for pressure in "//&
                 "boundary cells is extrapolated, rather than using PCM "//&
                 "in these cells. If true, the same order polynomial is "//&
                 "used as is used for the interior cells.", default=.true.)
 
  if (cs%tides) then
    cs%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
        time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=us%Z_to_m)
  endif
 
  cs%GFS_scale = 1.0
  if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
 
  call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale)
 
end subroutine pressureforce_afv_init
 
!> Deallocates the finite volume pressure gradient control structure
subroutine pressureforce_afv_end(CS)
  type(pressureforce_afv_cs), pointer :: cs !< Finite volume pressure control structure that
                                            !! will be deallocated in this subroutine.
  if (associated(cs)) deallocate(cs)
end subroutine pressureforce_afv_end
 
!> \namespace mom_pressureforce_afv
!!
!! Provides the Boussinesq and non-Boussinesq forms of horizontal accelerations
!! due to pressure gradients using a 2nd-order analytically vertically integrated
!! finite volume form, as described by Adcroft et al., 2008.
!!
!! This form eliminates the thermobaric instabilities that had been a problem with
!! previous forms of the pressure gradient force calculation, as described by
!! Hallberg, 2005.
!!
!! Adcroft, A., R. Hallberg, and M. Harrison, 2008: A finite volume discretization
!! of the pressure gradient force using analytic integration. Ocean Modelling, 22,
!! 106-113. http://doi.org/10.1016/j.ocemod.2008.02.001
!!
!! Hallberg, 2005: A thermobaric instability of Lagrangian vertical coordinate
!! ocean models. Ocean Modelling, 8, 279-300.
!! http://dx.doi.org/10.1016/j.ocemod.2004.01.001
 
end module mom_pressureforce_afv