mom_pressureforce_afv Module Reference

Detailed Description

Analytically integrated finite volume pressure gradient.

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

Data Types
type	pressureforce_afv_cs
	Finite volume pressure gradient control structure. More...

Functions/Subroutines
subroutine, public	pressureforce_afv (h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
	Thin interface between the model and the Boussinesq and non-Boussinesq pressure force routines. More...
subroutine, public	pressureforce_afv_nonbouss (h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
	Non-Boussinesq analytically-integrated finite volume form of pressure gradient. More...
subroutine, public	pressureforce_afv_bouss (h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
	Boussinesq analytically-integrated finite volume form of pressure gradient. More...
subroutine, public	pressureforce_afv_init (Time, G, GV, US, param_file, diag, CS, tides_CSp)
	Initializes the finite volume pressure gradient control structure. More...
subroutine, public	pressureforce_afv_end (CS)
	Deallocates the finite volume pressure gradient control structure. More...

Function/Subroutine Documentation

pressureforce_afv()

subroutine, public mom_pressureforce_afv::pressureforce_afv	(	real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(thermo_var_ptrs), intent(inout)	tv,
		real, dimension(szib_(g),szj_(g),szk_(g)), intent(out)	PFu,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)	PFv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(pressureforce_afv_cs), pointer	CS,
		type(ale_cs), pointer	ALE_CSp,
		real, dimension(:,:), optional, pointer	p_atm,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional	pbce,
		real, dimension(szi_(g),szj_(g)), intent(out), optional	eta
	)

Thin interface between the model and the Boussinesq and non-Boussinesq pressure force routines.

Parameters

[in]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thickness [H ~> m or kg m-2]
[in,out]	tv	Thermodynamic variables
[out]	pfu	Zonal acceleration [L T-2 ~> m s-2]
[out]	pfv	Meridional acceleration [L T-2 ~> m s-2]
	cs	Finite volume PGF control structure
	ale_csp	ALE control structure
	p_atm	The pressure at the ice-ocean or atmosphere-ocean interface [Pa].
[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].
[out]	eta	The bottom mass used to calculate PFu and PFv [H ~> m or kg m-2], with any tidal contributions or compressibility compensation.

Definition at line 67 of file MOM_PressureForce_analytic_FV.F90.

  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 [L T-2 ~> m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: PFv !< Meridional acceleration [L T-2 ~> 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
 

References pressureforce_afv_bouss(), and pressureforce_afv_nonbouss().

Here is the call graph for this function:

pressureforce_afv_bouss()

subroutine, public mom_pressureforce_afv::pressureforce_afv_bouss	(	real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		real, dimension(szib_(g),szj_(g),szk_(g)), intent(out)	PFu,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)	PFv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(pressureforce_afv_cs), pointer	CS,
		type(ale_cs), pointer	ALE_CSp,
		real, dimension(:,:), optional, pointer	p_atm,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional	pbce,
		real, dimension(szi_(g),szj_(g)), intent(out), optional	eta
	)

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).

Parameters

[in]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thickness [H ~> m]
[in]	tv	Thermodynamic variables
[out]	pfu	Zonal acceleration [L T-2 ~> m s-2]
[out]	pfv	Meridional acceleration [L T-2 ~> m s-2]
	cs	Finite volume PGF control structure
	ale_csp	ALE control structure
	p_atm	The pressure at the ice-ocean or atmosphere-ocean interface [Pa].
[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].
[out]	eta	The bottom mass used to calculate PFu and PFv [H ~> m or kg m-2], with any tidal contributions or compressibility compensation.

Definition at line 446 of file MOM_PressureForce_analytic_FV.F90.

  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 [L T-2 ~> m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> 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 [L2 T-2 ~> m2 s-2].
  real, dimension(SZI_(G)) :: &
    Rho_cv_BL   !   The coordinate potential density in the deepest variable
                ! density near-surface layer [R ~> kg m-3].
  real, dimension(SZI_(G),SZJ_(G)) :: &
    dz_geo, &   ! The change in geopotential thickness through a layer times some dimensional
                ! rescaling factors [kg m-1 R-1 s-2 ~> 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 [R ~> 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_mks_z      ! A scaled version of g_Earth [m2 Z-1 s-2 ~> m s-2].
  real :: g_Earth_z_geo      ! Another scaled version of g_Earth [R m5 kg-1 Z-1 s-2 ~> m s-2].
  real :: I_Rho0             ! 1/Rho0 times unit scaling factors [L2 m kg-1 s2 T-2 ~> m3 kg-1].
  real :: G_Rho0             ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
  real :: Rho_ref            ! The reference density [R ~> kg m-3].
  real :: Rho_ref_mks        ! The reference density in mks units [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 = us%m_s_to_L_T**2 / (us%R_to_kg_m3*gv%Rho0)
  g_earth_mks_z = us%L_T_to_m_s**2 * gv%g_Earth
  g_earth_z_geo = us%R_to_kg_m3*us%L_T_to_m_s**2 * gv%g_Earth
  g_rho0 = gv%g_Earth / gv%Rho0
  rho_ref_mks = cs%Rho0
  rho_ref = rho_ref_mks*us%kg_m3_to_R
 
  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, scale=us%kg_m3_to_R)
 
        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, scale=us%kg_m3_to_R)
        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, scale=us%kg_m3_to_R)
        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_geo)*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_geo)*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_mks, cs%Rho0, g_earth_mks_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_mks, cs%Rho0, g_earth_mks_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_mks, cs%Rho0, g_earth_mks_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_geo(i,j) = g_earth_z_geo * gv%H_to_Z*h(i,j,k)
        dpa(i,j) = (gv%Rlay(k) - rho_ref) * dz_geo(i,j)
        intz_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * dz_geo(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_geo(i,j) + dz_geo(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_geo(i,j) + dz_geo(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)
 

References mom_tidal_forcing::calc_tidal_forcing(), mom_error_handler::mom_error(), mom_ale::pressure_gradient_plm(), mom_ale::pressure_gradient_ppm(), and mom_pressureforce_mont::set_pbce_bouss().

Referenced by pressureforce_afv().

Here is the call graph for this function:

Here is the caller graph for this function:

pressureforce_afv_end()

subroutine, public mom_pressureforce_afv::pressureforce_afv_end

(

type(pressureforce_afv_cs), pointer

CS

)

Deallocates the finite volume pressure gradient control structure.

Parameters

cs	Finite volume pressure control structure that will be deallocated in this subroutine.

Definition at line 864 of file MOM_PressureForce_analytic_FV.F90.

  type(PressureForce_AFV_CS), pointer :: CS !< Finite volume pressure control structure that
                                            !! will be deallocated in this subroutine.
  if (associated(cs)) deallocate(cs)

pressureforce_afv_init()

subroutine, public mom_pressureforce_afv::pressureforce_afv_init	(	type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(pressureforce_afv_cs), pointer	CS,
		type(tidal_forcing_cs), optional, pointer	tides_CSp
	)

Initializes the finite volume pressure gradient control structure.

Parameters

[in]	time	Current model time
[in]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	param_file	Parameter file handles
[in,out]	diag	Diagnostics control structure
	cs	Finite volume PGF control structure
	tides_csp	Tides control structure

Definition at line 791 of file MOM_PressureForce_analytic_FV.F90.

  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)
 

References mom_error_handler::mom_error().

Here is the call graph for this function:

pressureforce_afv_nonbouss()

subroutine, public mom_pressureforce_afv::pressureforce_afv_nonbouss	(	real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		real, dimension(szib_(g),szj_(g),szk_(g)), intent(out)	PFu,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)	PFv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(pressureforce_afv_cs), pointer	CS,
		type(ale_cs), pointer	ALE_CSp,
		real, dimension(:,:), optional, pointer	p_atm,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(out), optional	pbce,
		real, dimension(szi_(g),szj_(g)), intent(out), optional	eta
	)

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).

Parameters

[in]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thickness [H ~> kg/m2]
[in]	tv	Thermodynamic variables
[out]	pfu	Zonal acceleration [L T-2 ~> m s-2]
[out]	pfv	Meridional acceleration [L T-2 ~> m s-2]
	cs	Finite volume PGF control structure
	ale_csp	ALE control structure
	p_atm	The pressure at the ice-ocean or atmosphere-ocean interface [Pa].
[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].
[out]	eta	The bottom mass used to calculate PFu and PFv [H ~> m or kg m-2], with any tidal contributions or compressibility compensation.

Definition at line 103 of file MOM_PressureForce_analytic_FV.F90.

  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 [L T-2 ~> m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> 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 [L2 T-2 ~> 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 [R ~> 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, scale=us%kg_m3_to_R)
        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/(us%R_to_kg_m3*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) * us%m_s_to_L_T**2 * &
            (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) * us%m_s_to_L_T**2 * &
          (p(i,j,1)*(1.0/(us%R_to_kg_m3*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))) * &
                   (us%m_s_to_L_T**2 * 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))) * &
                   (us%m_s_to_L_T**2 * 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)
 

References mom_tidal_forcing::calc_tidal_forcing(), mom_error_handler::mom_error(), mom_ale::pressure_gradient_plm(), mom_ale::pressure_gradient_ppm(), and mom_pressureforce_mont::set_pbce_nonbouss().

Referenced by pressureforce_afv().

Here is the call graph for this function:

Here is the caller graph for this function: