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, 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, CS, 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, 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, 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(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]	h	Layer thickness (m or kg/m2)
[in,out]	tv	Thermodynamic variables
[out]	pfu	Zonal acceleration (m/s2)
[out]	pfv	Meridional acceleration (m/s2)
	cs	Finite volume PGF control structure
	ale_csp	ALE control structure
	p_atm	The pressure at the ice-ocean or atmosphere-ocean interface in Pa.
[out]	pbce	The baroclinic pressure anomaly in each layer due to eta anomalies, in m2 s-2 H-1.
[out]	eta	The bottom mass used to calculate PFu and PFv, in H, with any tidal contributions or compressibility compensation.

Definition at line 51 of file MOM_PressureForce_analytic_FV.F90.

References pressureforce_afv_bouss(), and pressureforce_afv_nonbouss().

   type(ocean_grid_type),                     intent(in)    :: g   !< Ocean grid structure
   type(verticalgrid_type),                   intent(in)    :: gv  !< Vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h   !< Layer thickness (m or kg/m2)
   type(thermo_var_ptrs),                     intent(inout) :: tv  !< Thermodynamic variables
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out)   :: pfu !< Zonal acceleration (m/s2)
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: pfv !< Meridional acceleration (m/s2)
   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 in Pa.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  optional, intent(out) :: pbce !< The baroclinic pressure
                                                            !! anomaly in each layer due to eta anomalies,
                                                            !! in m2 s-2 H-1.
   real, dimension(SZI_(G),SZJ_(G)),          optional, intent(out) :: eta !< The bottom mass used to
                                                            !! calculate PFu and PFv, in H, with any tidal
                                                            !! contributions or compressibility compensation.
 
   if (gv%Boussinesq) then
     call pressureforce_afv_bouss(h, tv, pfu, pfv, g, gv, cs, ale_csp, p_atm, pbce, eta)
   else
     call pressureforce_afv_nonbouss(h, tv, pfu, pfv, g, gv, cs, p_atm, pbce, eta)
   endif
 

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(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 ie to ie, je to je range before this subroutine is called: h[ie+1] and h[je+1] and (if tveqn_of_state is set) T[ie+1], S[ie+1], T[je+1], and S[je+1].

Parameters

[in]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	h	Layer thickness (kg/m2)
[in]	tv	Thermodynamic variables
[out]	pfu	Zonal acceleration (m/s2)
[out]	pfv	Meridional acceleration (m/s2)
	cs	Finite volume PGF control structure
	ale_csp	ALE control structure
	p_atm	The pressure at the ice-ocean or atmosphere-ocean interface in Pa.
[out]	pbce	The baroclinic pressure anomaly in each layer due to eta anomalies, in m2 s-2 H-1.
[out]	eta	The bottom mass used to calculate PFu and PFv, in H, with any tidal contributions or compressibility compensation.

Definition at line 414 of file MOM_PressureForce_analytic_FV.F90.

References mom_tidal_forcing::calc_tidal_forcing(), mom_error_handler::mom_error(), regrid_defs::pressure_reconstruction_plm, regrid_defs::pressure_reconstruction_ppm, and mom_pressureforce_mont::set_pbce_bouss().

Referenced by pressureforce_afv().

   type(ocean_grid_type),                     intent(in)  :: g   !< Ocean grid structure
   type(verticalgrid_type),                   intent(in)  :: gv  !< Vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)  :: h   !< Layer thickness (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/s2)
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration (m/s2)
   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 in Pa.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  optional, intent(out) :: pbce !< The baroclinic pressure
                                                          !! anomaly in each layer due to eta anomalies,
                                                          !! in m2 s-2 H-1.
   real, dimension(SZI_(G),SZJ_(G)),          optional, intent(out) :: eta !< The bottom mass used to
                                                          !! calculate PFu and PFv, in H, with any tidal
                                                          !! contributions or compressibility compensation.
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in 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, in m.
     dm          ! The barotropic adjustment to the Montgomery potential to
                 ! account for a reduced gravity model, in m2 s-2.
   real, dimension(SZI_(G)) :: &
     rho_cv_bl   !   The coordinate potential density in the deepest variable
                 ! density near-surface layer, in kg m-3.
   real, dimension(SZDI_(G%Block(1)),SZDJ_(G%Block(1))) :: &  ! on block indices
     dz_bk, &     ! The change in geopotential thickness through a layer, m2 s-2.
     pa_bk, &     ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the
                  ! the interface atop a layer, in Pa.
     dpa_bk, &    ! The change in pressure anomaly between the top and bottom
                  ! of a layer, in Pa.
     intz_dpa_bk  ! The vertical integral in depth of the pressure anomaly less
                  ! the pressure anomaly at the top of the layer, in H Pa (m Pa).
   real, dimension(SZDIB_(G%Block(1)),SZDJ_(G%Block(1))) :: & ! on block indices
     intx_pa_bk, & ! The zonal integral of the pressure anomaly along the interface
                   ! atop a layer, divided by the grid spacing, in Pa.
     intx_dpa_bk   ! The change in intx_pa through a layer, in Pa.
   real, dimension(SZDI_(G%Block(1)),SZDJB_(G%Block(1))) :: & ! on block indices
     inty_pa_bk, & ! The meridional integral of the pressure anomaly along the
                   ! interface atop a layer, divided by the grid spacing, in Pa.
     inty_dpa_bk   ! The change in inty_pa through a layer, in 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, in C.
     s_tmp       ! Temporary array of salinities where layers that are lighter
                 ! than the mixed layer have the mixed layer's properties, in psu.
   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, in kg m-3.
   real :: p_ref(szi_(g))     !   The pressure used to calculate the coordinate
                              ! density, in Pa (usually 2e7 Pa = 2000 dbar).
   real :: p0(szi_(g)) ! An array of zeros to use for pressure in Pa.
   real :: h_neglect          ! A thickness that is so small it is usually lost
                              ! in roundoff and can be neglected, in m.
   real :: i_rho0             ! 1/Rho0.
   real :: g_rho0             ! G_Earth / Rho0 in m4 s-2 kg-1.
   real :: rho_ref            ! The reference density in kg m-3.
   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 :: is_bk, ie_bk, js_bk, je_bk, isq_bk, ieq_bk, jsq_bk, jeq_bk
   integer :: ioff_bk, joff_bk
   integer :: i, j, k, n, ib, jb
   integer :: prscheme
 
   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: 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 = usepressurereconstruction(ale_csp) .and. use_eos
 
   prscheme = pressurereconstructionscheme(ale_csp)
   h_neglect = gv%H_subroundoff
   i_rho0 = 1.0/gv%Rho0
   g_rho0 = gv%g_Earth/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(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,GV,h)
     do j=jsq,jeq+1
       do i=isq,ieq+1
         e(i,j,1) = -1.0*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_m
       enddo ; enddo
     enddo
     call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp)
   endif
 
 !    Here layer interface heights, e, are calculated.
 !$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,GV,h,CS,e_tidal)
   if (cs%tides) then
 !$OMP do
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       e(i,j,nz+1) = -1.0*g%bathyT(i,j) - e_tidal(i,j)
     enddo ; enddo
   else
 !$OM do
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       e(i,j,nz+1) = -1.0*g%bathyT(i,j)
     enddo ; enddo
   endif
 !$OMP do
   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_m
   enddo ; enddo ; enddo
 !$OMP end parallel
 
 
   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(none) shared(Isq,Ieq,Jsq,Jeq,nkmb,nz,GV,tv_tmp,tv,p_ref) &
 !$OMP                          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
 
 !$OMP parallel default(none) shared(Jsq,Jeq,Isq,Ieq,tv_tmp,p_atm,rho_in_situ,tv, &
 !$OMP                               p0,dM,CS,G_Rho0,e,use_p_atm,use_EOS,GV,    &
 !$OMP                               rho_ref,js,je,is,ie)
   if (cs%GFS_scale < 1.0) then
     ! Adjust the Montgomery potential to make this a reduced gravity model.
     if (use_eos) then
 !$OMP do
       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 do
       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
 !$OMP end parallel
 
 ! 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 ( prscheme == pressure_reconstruction_plm ) then
       call pressure_gradient_plm (ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h);
     elseif ( prscheme == pressure_reconstruction_ppm ) then
       call pressure_gradient_ppm (ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h);
     endif
   endif
 
 !$OMP parallel do default(none) shared(use_p_atm,rho_ref,G,GV,e,     &
 !$OMP                                  p_atm,nz,use_EOS,use_ALE,PRScheme,T_t,T_b,S_t, &
 !$OMP                                  S_b,CS,tv,tv_tmp,h,PFu,I_Rho0,h_neglect,PFv,dM)&
 !$OMP                          private(is_bk,ie_bk,js_bk,je_bk,Isq_bk,Ieq_bk,Jsq_bk,  &
 !$OMP                                  Jeq_bk,ioff_bk,joff_bk,pa_bk,  &
 !$OMP                                  intx_pa_bk,inty_pa_bk,dpa_bk,intz_dpa_bk,      &
 !$OMP                                  intx_dpa_bk,inty_dpa_bk,dz_bk,i,j)
   do n = 1, g%nblocks
     is_bk=g%Block(n)%isc      ; ie_bk=g%Block(n)%iec
     js_bk=g%Block(n)%jsc      ; je_bk=g%Block(n)%jec
     isq_bk=g%Block(n)%IscB    ; ieq_bk=g%Block(n)%IecB
     jsq_bk=g%Block(n)%JscB    ; jeq_bk=g%Block(n)%JecB
     ioff_bk = g%Block(n)%idg_offset - g%HI%idg_offset
     joff_bk = g%Block(n)%jdg_offset - g%HI%jdg_offset
 
     ! 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
       do jb=jsq_bk,jeq_bk+1 ; do ib=isq_bk,ieq_bk+1
         i = ib+ioff_bk ; j = jb+joff_bk
         pa_bk(ib,jb) = (rho_ref*gv%g_Earth)*e(i,j,1) + p_atm(i,j)
       enddo ; enddo
     else
       do jb=jsq_bk,jeq_bk+1 ; do ib=isq_bk,ieq_bk+1
         i = ib+ioff_bk ; j = jb+joff_bk
         pa_bk(ib,jb) = (rho_ref*gv%g_Earth)*e(i,j,1)
       enddo ; enddo
     endif
     do jb=js_bk,je_bk ; do ib=isq_bk,ieq_bk
       intx_pa_bk(ib,jb) = 0.5*(pa_bk(ib,jb) + pa_bk(ib+1,jb))
     enddo ; enddo
     do jb=jsq_bk,jeq_bk ; do ib=is_bk,ie_bk
       inty_pa_bk(ib,jb) = 0.5*(pa_bk(ib,jb) + pa_bk(ib,jb+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 ( prscheme == pressure_reconstruction_plm ) 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, gv%g_Earth,    &
                       gv%H_subroundoff, g%bathyT, g%HI, g%Block(n), &
                       tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, &
                       usemasswghtinterp = cs%useMassWghtInterp)
           elseif ( prscheme == pressure_reconstruction_ppm ) 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, gv%g_Earth, &
                       g%HI, g%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk,    &
                       intx_dpa_bk, inty_dpa_bk)
           endif
         else
           call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), &
                     e(:,:,k), e(:,:,k+1),             &
                     rho_ref, cs%Rho0, gv%g_Earth, g%HI, g%Block(n), tv%eqn_of_state, &
                     dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk )
         endif
         do jb=jsq_bk,jeq_bk+1 ; do ib=isq_bk,ieq_bk+1
           intz_dpa_bk(ib,jb) = intz_dpa_bk(ib,jb)*gv%m_to_H
         enddo ; enddo
       else
         do jb=jsq_bk,jeq_bk+1 ; do ib=isq_bk,ieq_bk+1
           i = ib+ioff_bk ; j = jb+joff_bk
           dz_bk(ib,jb) = gv%g_Earth*gv%H_to_m*h(i,j,k)
           dpa_bk(ib,jb) = (gv%Rlay(k) - rho_ref)*dz_bk(ib,jb)
           intz_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref)*dz_bk(ib,jb)*h(i,j,k)
         enddo ; enddo
         do jb=js_bk,je_bk ; do ib=isq_bk,ieq_bk
           intx_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib+1,jb))
         enddo ; enddo
         do jb=jsq_bk,jeq_bk ; do ib=is_bk,ie_bk
           inty_dpa_bk(ib,jb) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib,jb+1))
         enddo ; enddo
       endif
 
       ! Compute pressure gradient in x direction
       do jb=js_bk,je_bk ; do ib=isq_bk,ieq_bk
         i = ib+ioff_bk ; j = jb+joff_bk
         pfu(i,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
                      (pa_bk(ib+1,jb)*h(i+1,j,k) + intz_dpa_bk(ib+1,jb))) + &
                      ((h(i+1,j,k) - h(i,j,k)) * intx_pa_bk(ib,jb) - &
                      (e(i+1,j,k+1) - e(i,j,k+1)) * intx_dpa_bk(ib,jb) * gv%m_to_H)) * &
                      ((2.0*i_rho0*g%IdxCu(i,j)) / &
                      ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
         intx_pa_bk(ib,jb) = intx_pa_bk(ib,jb) + intx_dpa_bk(ib,jb)
       enddo ; enddo
       ! Compute pressure gradient in y direction
       do jb=jsq_bk,jeq_bk ; do ib=is_bk,ie_bk
         i = ib+ioff_bk ; j = jb+joff_bk
         pfv(i,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
                      (pa_bk(ib,jb+1)*h(i,j+1,k) + intz_dpa_bk(ib,jb+1))) + &
                      ((h(i,j+1,k) - h(i,j,k)) * inty_pa_bk(ib,jb) - &
                      (e(i,j+1,k+1) - e(i,j,k+1)) * inty_dpa_bk(ib,jb) * gv%m_to_H)) * &
                      ((2.0*i_rho0*g%IdyCv(i,j)) / &
                      ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
         inty_pa_bk(ib,jb) = inty_pa_bk(ib,jb) + inty_dpa_bk(ib,jb)
       enddo ; enddo
       do jb=jsq_bk,jeq_bk+1 ; do ib=isq_bk,ieq_bk+1
         pa_bk(ib,jb) = pa_bk(ib,jb) + dpa_bk(ib,jb)
       enddo ; enddo
     enddo
 
     if (cs%GFS_scale < 1.0) then
       do k=1,nz
         do j=js_bk+joff_bk,je_bk+joff_bk ; do i=isq_bk+ioff_bk,ieq_bk+ioff_bk
           pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
         enddo ; enddo
         do j=jsq_bk+joff_bk,jeq_bk+joff_bk ; do i=is_bk+ioff_bk,ie_bk+ioff_bk
           pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
         enddo ; enddo
       enddo
     endif
   enddo
 
   if (present(pbce)) then
     call set_pbce_bouss(e, tv_tmp, g, gv, gv%g_Earth, 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.
 !$OM parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e,e_tidal)
       do j=jsq,jeq+1 ; do i=isq,ieq+1
         eta(i,j) = e(i,j,1)*gv%m_to_H + e_tidal(i,j)*gv%m_to_H
       enddo ; enddo
     else
 !$OM parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e)
       do j=jsq,jeq+1 ; do i=isq,ieq+1
         eta(i,j) = e(i,j,1)*gv%m_to_H
       enddo ; enddo
     endif
   endif
 
   if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
 

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.

Definition at line 821 of file MOM_PressureForce_analytic_FV.F90.

821 type(pressureforce_afv_cs), pointer :: cs

822 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(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]	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 771 of file MOM_PressureForce_analytic_FV.F90.

References mom_error_handler::mom_error().

Referenced by mom_pressureforce::pressureforce_init().

   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(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.
 
   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 \n"//&
                  "calculate accelerations and the mass for conservation \n"//&
                  "properties, or with BOUSSINSEQ false to convert some \n"//&
                  "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, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", cs%useMassWghtInterp, &
                  "If true, use mass weighting when interpolation T/S for\n"//&
                  "top/bottom integrals in AFV pressure gradient calculation.", default=.false.)
 
   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')
   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)
 

Here is the call graph for this function:

Here is the caller 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(pressureforce_afv_cs), pointer	CS,
		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 ie to ie, je to je range before this subroutine is called: h[ie+1] and h[je+1] and (if tveqn_of_state is set) T[ie+1], S[ie+1], T[je+1], and S[je+1].

Parameters

[in]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	h	Layer thickness (kg/m2)
[in]	tv	Thermodynamic variables
[out]	pfu	Zonal acceleration (m/s2)
[out]	pfv	Meridional acceleration (m/s2)
	cs	Finite volume PGF control structure
	p_atm	The pressure at the ice-ocean or atmosphere-ocean interface in Pa.
[out]	pbce	The baroclinic pressure anomaly in each layer due to eta anomalies, in m2 s-2 H-1.
[out]	eta	The bottom mass used to calculate PFu and PFv, in H, with any tidal contributions or compressibility compensation.

Definition at line 87 of file MOM_PressureForce_analytic_FV.F90.

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

Referenced by pressureforce_afv().

   type(ocean_grid_type),                     intent(in)    :: g   !< Ocean grid structure
   type(verticalgrid_type),                   intent(in)    :: gv  !< Vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h   !< Layer thickness (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/s2)
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: pfv !< Meridional acceleration (m/s2)
   type(pressureforce_afv_cs),                pointer       :: cs  !< Finite volume PGF control structure
   real, dimension(:,:),                      optional, pointer :: p_atm !< The pressure at the ice-ocean
                                                            !! or atmosphere-ocean interface in Pa.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  optional, intent(out) :: pbce !< The baroclinic pressure
                                                            !! anomaly in each layer due to eta anomalies,
                                                            !! in m2 s-2 H-1.
   real, dimension(SZI_(G),SZJ_(G)),          optional, intent(out) :: eta !< The bottom mass used to
                                                            !! calculate PFu and PFv, in H, with any tidal
                                                            !! contributions or compressibility compensation.
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure in 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, in C.
     s_tmp       ! Temporary array of salinities where layers that are lighter
                 ! than the mixed layer have the mixed layer's properties, in psu.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: &
     dza, &      ! The change in geopotential anomaly between the top and bottom
                 ! of a layer, in m2 s-2.
     intp_dza    ! The vertical integral in depth of the pressure anomaly less
                 ! the pressure anomaly at the top of the layer, in Pa m2 s-2.
   real, dimension(SZI_(G),SZJ_(G))  :: &
     dp, &       ! The (positive) change in pressure across a layer, in Pa.
     ssh, &      ! The sea surface height anomaly, in m.
     e_tidal, &  ! The bottom geopotential anomaly due to tidal forces from
                 ! astronomical sources and self-attraction and loading, in m.
     dm, &       ! The barotropic adjustment to the Montgomery potential to
                 ! account for a reduced gravity model, in m2 s-2.
     za          ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the
                 ! interface atop a layer, in m2 s-2.
   real, dimension(SZDI_(G%Block(1)),SZDJ_(G%Block(1))) :: & ! on block indices
     dp_bk, &    ! The (positive) change in pressure across a layer, in Pa.
     za_bk       ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the
                 ! interface atop a layer, in m2 s-2.
 
   real, dimension(SZI_(G)) :: rho_cv_bl !  The coordinate potential density in the deepest variable
                 ! density near-surface layer, in kg m-3.
   real, dimension(SZDIB_(G%Block(1)),SZDJ_(G%Block(1))) :: & ! on block indices
     intx_za_bk ! 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, in m2 s-2.
   real, dimension(SZDI_(G%Block(1)),SZDJB_(G%Block(1))) :: & ! on block indices
     inty_za_bk ! 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, in m2 s-2.
   real :: p_ref(szi_(g))     !   The pressure used to calculate the coordinate
                              ! density, in 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, in Pa.
   real :: alpha_anom         ! The in-situ specific volume, averaged over a
                              ! layer, less alpha_ref, in m3 kg-1.
   logical :: use_p_atm       ! If true, use the atmospheric pressure.
   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, in m3 kg-1, that is used
                     ! to reduce the impact of truncation errors.
   real :: rho_in_situ(szi_(g)) ! The in situ density, in 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/ms2)
   real :: i_gearth
   real, parameter :: c1_6 = 1.0/6.0
   integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
   integer :: is_bk, ie_bk, js_bk, je_bk, isq_bk, ieq_bk, jsq_bk, jeq_bk
   integer :: i, j, k, n, ib, jb, ioff_bk, joff_bk
 
   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
 
   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)
 
   if (.not.associated(cs)) call mom_error(fatal, &
        "MOM_PressureForce: Module must be initialized before it is used.")
 
   dp_neglect = gv%H_to_Pa * gv%H_subroundoff
   alpha_ref = 1.0/cs%Rho0
 
 !$OMP parallel default(none) shared(Isq,Ieq,Jsq,Jeq,nz,use_p_atm,p,p_atm,GV,h)
   if (use_p_atm) then
 !$OMP do
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       p(i,j,1) = p_atm(i,j)
     enddo ; enddo
   else
 !$OMP do
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       p(i,j,1) = 0.0 ! or oneatm
     enddo ; enddo
   endif
 !$OMP do
   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
 !$OMP end parallel
 
   i_gearth = 1.0 / gv%g_Earth
 
   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(none) shared(Isq,Ieq,Jsq,Jeq,nz,nkmb,tv_tmp,tv,p_ref,GV) &
 !$OMP                          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
 
 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,is,ie,js,je,tv_tmp,alpha_ref, &
 !$OMP                                  p,h,G,GV,tv,dza,intp_dza,intx_dza,inty_dza,use_EOS) &
 !$OMP                          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
       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))
     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(none) shared(Isq,Ieq,Jsq,Jeq,nz,za,alpha_ref,p,G,GV,dza)
   do j=jsq,jeq+1
     do i=isq,ieq+1
       za(i,j) = alpha_ref*p(i,j,nz+1) - gv%g_Earth*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(none) shared(Isq,Ieq,Jsq,Jeq,SSH,za,alpha_ref,p,I_gEarth)
     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)
 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,za,G,GV,e_tidal)
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       za(i,j) = za(i,j) - gv%g_Earth*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(none) shared(Isq,Ieq,Jsq,Jeq,tv_tmp,p,tv,dM,CS,alpha_ref,za) &
 !$OMP                          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(none) shared(Isq,Ieq,Jsq,Jeq,dM,CS,p,GV,alpha_ref,za)
       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(none) shared(nz,za,G,GV,dza,intx_dza,h,PFu, &
 !$OMP                                  intp_dza,p,dp_neglect,inty_dza,PFv,CS,dM) &
 !$OMP                          private(is_bk,ie_bk,js_bk,je_bk,Isq_bk,Ieq_bk,Jsq_bk, &
 !$OMP                                  Jeq_bk,ioff_bk,joff_bk,i,j,za_bk,intx_za_bk,  &
 !$OMP                                  inty_za_bk,dp_bk)
   do n = 1, g%nblocks
     is_bk=g%block(n)%isc      ; ie_bk=g%block(n)%iec
     js_bk=g%block(n)%jsc      ; je_bk=g%block(n)%jec
     isq_bk=g%block(n)%IscB    ; ieq_bk=g%block(n)%IecB
     jsq_bk=g%block(n)%JscB    ; jeq_bk=g%block(n)%JecB
     ioff_bk = g%Block(n)%idg_offset - g%HI%idg_offset
     joff_bk = g%Block(n)%jdg_offset - g%HI%jdg_offset
     do jb=jsq_bk,jeq_bk+1 ; do ib=isq_bk,ieq_bk+1
       i = ib+ioff_bk ; j = jb+joff_bk
       za_bk(ib,jb) = za(i,j)
     enddo ; enddo
     do jb=js_bk,je_bk ; do ib=isq_bk,ieq_bk
       i = ib+ioff_bk ; j = jb+joff_bk
       intx_za_bk(ib,jb) = 0.5*(za_bk(ib,jb) + za_bk(ib+1,jb))
     enddo ; enddo
     do jb=jsq_bk,jeq_bk ; do ib=is_bk,ie_bk
       i = ib+ioff_bk ; j = jb+joff_bk
       inty_za_bk(ib,jb) = 0.5*(za_bk(ib,jb) + za_bk(ib,jb+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.
       do jb=jsq_bk,jeq_bk+1 ; do ib=isq_bk,ieq_bk+1
         i = ib+ioff_bk ; j = jb+joff_bk
         dp_bk(ib,jb) = gv%H_to_Pa*h(i,j,k)
         za_bk(ib,jb) = za_bk(ib,jb) - dza(i,j,k)
       enddo ; enddo
       do jb=js_bk,je_bk ; do ib=isq_bk,ieq_bk
         i = ib+ioff_bk ; j = jb+joff_bk
         intx_za_bk(ib,jb) = intx_za_bk(ib,jb) - intx_dza(i,j,k)
         pfu(i,j,k) = (((za_bk(ib,jb)*dp_bk(ib,jb) + intp_dza(i,j,k)) - &
                      (za_bk(ib+1,jb)*dp_bk(ib+1,jb) + intp_dza(i+1,j,k))) + &
                      ((dp_bk(ib+1,jb) - dp_bk(ib,jb)) * intx_za_bk(ib,jb) - &
                      (p(i+1,j,k) - p(i,j,k)) * intx_dza(i,j,k))) * &
                      (2.0*g%IdxCu(i,j) / ((dp_bk(ib,jb) + dp_bk(ib+1,jb)) + &
                      dp_neglect))
       enddo ; enddo
       do jb=jsq_bk,jeq_bk ; do ib=is_bk,ie_bk
         i = ib+ioff_bk ; j = jb+joff_bk
         inty_za_bk(ib,jb) = inty_za_bk(ib,jb) - inty_dza(i,j,k)
         pfv(i,j,k) = (((za_bk(ib,jb)*dp_bk(ib,jb) + intp_dza(i,j,k)) - &
                      (za_bk(ib,jb+1)*dp_bk(ib,jb+1) + intp_dza(i,j+1,k))) + &
                      ((dp_bk(ib,jb+1) - dp_bk(ib,jb)) * inty_za_bk(ib,jb) - &
                      (p(i,j+1,k) - p(i,j,k)) * inty_dza(i,j,k))) * &
                      (2.0*g%IdyCv(i,j) / ((dp_bk(ib,jb) + dp_bk(ib,jb+1)) + &
                      dp_neglect))
       enddo ; enddo
 
       if (cs%GFS_scale < 1.0) then
         ! Adjust the Montgomery potential to make this a reduced gravity model.
         do j=js_bk+joff_bk,je_bk+joff_bk ; do i=isq_bk+ioff_bk,ieq_bk+ioff_bk
           pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
         enddo ; enddo
         do j=jsq_bk+joff_bk,jeq_bk+joff_bk ; do i=is_bk+ioff_bk,ie_bk+ioff_bk
           pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
         enddo ; enddo
       endif
     enddo
   enddo
 
   if (present(pbce)) then
     call set_pbce_nonbouss(p, tv_tmp, g, gv, gv%g_Earth, 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(none) shared(Isq,Ieq,Jsq,Jeq,nz,eta,p,p_atm,Pa_to_H)
       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(none) shared(Isq,Ieq,Jsq,Jeq,nz,eta,p,Pa_to_H)
       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)
 

Here is the call graph for this function:

Here is the caller graph for this function:

Detailed Description

Data Types

Functions/Subroutines

Function/Subroutine Documentation

◆ pressureforce_afv()

◆ pressureforce_afv_bouss()

◆ pressureforce_afv_end()

◆ pressureforce_afv_init()

◆ pressureforce_afv_nonbouss()