MOM6
Data Types | Functions/Subroutines
mom_pressureforce_fv Module Reference

Detailed Description

Finite volume pressure gradient (integrated by quadrature or analytically)

Provides the Boussinesq and non-Boussinesq forms of horizontal accelerations due to pressure gradients using a vertically integrated finite volume form, as described by Adcroft et al., 2008. Integration in the vertical is made either by quadrature or analytically.

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_fv_cs
 Finite volume pressure gradient control structure. More...
 

Functions/Subroutines

subroutine, public pressureforce_fv_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_fv_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_fv_init (Time, G, GV, US, param_file, diag, CS, tides_CSp)
 Initializes the finite volume pressure gradient control structure. More...
 
subroutine, public pressureforce_fv_end (CS)
 Deallocates the finite volume pressure gradient control structure. More...
 

Function/Subroutine Documentation

◆ pressureforce_fv_bouss()

subroutine, public mom_pressureforce_fv::pressureforce_fv_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_fv_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]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thickness [H ~> m]
[in]tvThermodynamic variables
[out]pfuZonal acceleration [L T-2 ~> m s-2]
[out]pfvMeridional acceleration [L T-2 ~> m s-2]
csFinite volume PGF control structure
ale_cspALE control structure
p_atmThe pressure at the ice-ocean or atmosphere-ocean interface [R L2 T-2 ~> Pa].
[out]pbceThe baroclinic pressure anomaly in each layer due to eta anomalies [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
[out]etaThe bottom mass used to calculate PFu and PFv [H ~> m or kg m-2], with any tidal contributions or compressibility compensation.

Definition at line 416 of file MOM_PressureForce_FV.F90.

417  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
418  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
419  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
420  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m]
421  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
422  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: PFu !< Zonal acceleration [L T-2 ~> m s-2]
423  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> m s-2]
424  type(PressureForce_FV_CS), pointer :: CS !< Finite volume PGF control structure
425  type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
426  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean
427  !! or atmosphere-ocean interface [R L2 T-2 ~> Pa].
428  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
429  !! anomaly in each layer due to eta anomalies
430  !! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
431  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
432  !! calculate PFu and PFv [H ~> m or kg m-2], with any
433  !! tidal contributions or compressibility compensation.
434  ! Local variables
435  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in depth units [Z ~> m].
436  real, dimension(SZI_(G),SZJ_(G)) :: &
437  e_tidal, & ! The bottom geopotential anomaly due to tidal forces from
438  ! astronomical sources and self-attraction and loading [Z ~> m].
439  dm ! The barotropic adjustment to the Montgomery potential to
440  ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
441  real, dimension(SZI_(G)) :: &
442  Rho_cv_BL ! The coordinate potential density in the deepest variable
443  ! density near-surface layer [R ~> kg m-3].
444  real, dimension(SZI_(G),SZJ_(G)) :: &
445  dz_geo, & ! The change in geopotential thickness through a layer [L2 T-2 ~> m2 s-2].
446  pa, & ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the
447  ! the interface atop a layer [R L2 T-2 ~> Pa].
448  dpa, & ! The change in pressure anomaly between the top and bottom
449  ! of a layer [R L2 T-2 ~> Pa].
450  intz_dpa ! The vertical integral in depth of the pressure anomaly less the
451  ! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa or kg m-2 Pa].
452  real, dimension(SZIB_(G),SZJ_(G)) :: &
453  intx_pa, & ! The zonal integral of the pressure anomaly along the interface
454  ! atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].
455  intx_dpa ! The change in intx_pa through a layer [R L2 T-2 ~> Pa].
456  real, dimension(SZI_(G),SZJB_(G)) :: &
457  inty_pa, & ! The meridional integral of the pressure anomaly along the
458  ! interface atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].
459  inty_dpa ! The change in inty_pa through a layer [R L2 T-2 ~> Pa].
460 
461  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
462  T_tmp, & ! Temporary array of temperatures where layers that are lighter
463  ! than the mixed layer have the mixed layer's properties [degC].
464  s_tmp ! Temporary array of salinities where layers that are lighter
465  ! than the mixed layer have the mixed layer's properties [ppt].
466  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
467  S_t, S_b, T_t, T_b ! Top and bottom edge values for linear reconstructions
468  ! of salinity and temperature within each layer.
469  real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3].
470  real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
471  ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
472  real :: p0(SZI_(G)) ! An array of zeros to use for pressure [R L2 T-2 ~> Pa].
473  real :: h_neglect ! A thickness that is so small it is usually lost
474  ! in roundoff and can be neglected [H ~> m].
475  real :: I_Rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1].
476  real :: G_Rho0 ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
477  real :: rho_ref ! The reference density [R ~> kg m-3].
478  real :: dz_neglect ! A minimal thickness [Z ~> m], like e.
479  logical :: use_p_atm ! If true, use the atmospheric pressure.
480  logical :: use_ALE ! If true, use an ALE pressure reconstruction.
481  logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
482  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
483  real :: Tl(5) ! copy and T in local stencil [degC]
484  real :: mn_T ! mean of T in local stencil [degC]
485  real :: mn_T2 ! mean of T**2 in local stencil [degC]
486  real :: hl(5) ! Copy of local stencil of H [H ~> m]
487  real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1]
488  real, parameter :: C1_6 = 1.0/6.0
489  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
490  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
491  integer :: i, j, k
492 
493  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
494  nkmb=gv%nk_rho_varies
495  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
496  eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)
497 
498  if (.not.associated(cs)) call mom_error(fatal, &
499  "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.")
500 
501  use_p_atm = .false.
502  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
503  use_eos = associated(tv%eqn_of_state)
504  do i=isq,ieq+1 ; p0(i) = 0.0 ; enddo
505  use_ale = .false.
506  if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
507 
508  if (cs%Stanley_T2_det_coeff>=0.) then
509  if (.not. associated(tv%varT)) call safe_alloc_ptr(tv%varT, g%isd, g%ied, g%jsd, g%jed, gv%ke)
510  do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1
511  ! Strictly speaking we should estimate the *horizontal* grid-scale variance
512  ! but neither of the following blocks make a rotation to the horizontal
513  ! and instead work along coordinate.
514 
515  ! This block calculates a simple |delta T| along coordinates and does
516  ! not allow vanishing layer thicknesses or layers tracking topography
517  !! SGS variance in i-direction [degC2]
518  !dTdi2 = ( ( G%mask2dCu(I ,j) * G%IdxCu(I ,j) * ( tv%T(i+1,j,k) - tv%T(i,j,k) ) &
519  ! + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( tv%T(i,j,k) - tv%T(i-1,j,k) ) &
520  ! ) * G%dxT(i,j) * 0.5 )**2
521  !! SGS variance in j-direction [degC2]
522  !dTdj2 = ( ( G%mask2dCv(i,J ) * G%IdyCv(i,J ) * ( tv%T(i,j+1,k) - tv%T(i,j,k) ) &
523  ! + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( tv%T(i,j,k) - tv%T(i,j-1,k) ) &
524  ! ) * G%dyT(i,j) * 0.5 )**2
525  !tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * 0.5 * ( dTdi2 + dTdj2 )
526 
527  ! This block does a thickness weighted variance calculation and helps control for
528  ! extreme gradients along layers which are vanished against topography. It is
529  ! still a poor approximation in the interior when coordinates are strongly tilted.
530  hl(1) = h(i,j,k) * g%mask2dT(i,j)
531  hl(2) = h(i-1,j,k) * g%mask2dCu(i-1,j)
532  hl(3) = h(i+1,j,k) * g%mask2dCu(i,j)
533  hl(4) = h(i,j-1,k) * g%mask2dCv(i,j-1)
534  hl(5) = h(i,j+1,k) * g%mask2dCv(i,j)
535  r_sm_h = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + gv%H_subroundoff )
536  ! Mean of T
537  tl(1) = tv%T(i,j,k) ; tl(2) = tv%T(i-1,j,k) ; tl(3) = tv%T(i+1,j,k)
538  tl(4) = tv%T(i,j-1,k) ; tl(5) = tv%T(i,j+1,k)
539  mn_t = ( hl(1)*tl(1) + ( ( hl(2)*tl(2) + hl(3)*tl(3) ) + ( hl(4)*tl(4) + hl(5)*tl(5) ) ) ) * r_sm_h
540  ! Adjust T vectors to have zero mean
541  tl(:) = tl(:) - mn_t ; mn_t = 0.
542  ! Variance of T
543  mn_t2 = ( hl(1)*tl(1)*tl(1) + ( ( hl(2)*tl(2)*tl(2) + hl(3)*tl(3)*tl(3) ) &
544  + ( hl(4)*tl(4)*tl(4) + hl(5)*tl(5)*tl(5) ) ) ) * r_sm_h
545  ! Variance should be positive but round-off can violate this. Calculating
546  ! variance directly would fix this but requires more operations.
547  tv%varT(i,j,k) = cs%Stanley_T2_det_coeff * max(0., mn_t2)
548  enddo ; enddo ; enddo
549  endif
550 
551  h_neglect = gv%H_subroundoff
552  dz_neglect = gv%H_subroundoff * gv%H_to_Z
553  i_rho0 = 1.0 / gv%Rho0
554  g_rho0 = gv%g_Earth / gv%Rho0
555  rho_ref = cs%Rho0
556 
557  if (cs%tides) then
558  ! Determine the surface height anomaly for calculating self attraction
559  ! and loading. This should really be based on bottom pressure anomalies,
560  ! but that is not yet implemented, and the current form is correct for
561  ! barotropic tides.
562  !$OMP parallel do default(shared)
563  do j=jsq,jeq+1
564  do i=isq,ieq+1
565  e(i,j,1) = -g%bathyT(i,j)
566  enddo
567  do k=1,nz ; do i=isq,ieq+1
568  e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
569  enddo ; enddo
570  enddo
571  call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
572  endif
573 
574 ! Here layer interface heights, e, are calculated.
575  if (cs%tides) then
576  !$OMP parallel do default(shared)
577  do j=jsq,jeq+1 ; do i=isq,ieq+1
578  e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
579  enddo ; enddo
580  else
581  !$OMP parallel do default(shared)
582  do j=jsq,jeq+1 ; do i=isq,ieq+1
583  e(i,j,nz+1) = -g%bathyT(i,j)
584  enddo ; enddo
585  endif
586  !$OMP parallel do default(shared)
587  do j=jsq,jeq+1; do k=nz,1,-1 ; do i=isq,ieq+1
588  e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
589  enddo ; enddo ; enddo
590 
591  if (use_eos) then
592 ! With a bulk mixed layer, replace the T & S of any layers that are
593 ! lighter than the the buffer layer with the properties of the buffer
594 ! layer. These layers will be massless anyway, and it avoids any
595 ! formal calculations with hydrostatically unstable profiles.
596 
597  if (nkmb>0) then
598  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
599  tv_tmp%eqn_of_state => tv%eqn_of_state
600 
601  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
602  !$OMP parallel do default(shared) private(Rho_cv_BL)
603  do j=jsq,jeq+1
604  do k=1,nkmb ; do i=isq,ieq+1
605  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
606  enddo ; enddo
607  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, rho_cv_bl(:), &
608  tv%eqn_of_state, eosdom)
609 
610  do k=nkmb+1,nz ; do i=isq,ieq+1
611  if (gv%Rlay(k) < rho_cv_bl(i)) then
612  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
613  else
614  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
615  endif
616  enddo ; enddo
617  enddo
618  else
619  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
620  tv_tmp%eqn_of_state => tv%eqn_of_state
621  endif
622  endif
623 
624  if (cs%GFS_scale < 1.0) then
625  ! Adjust the Montgomery potential to make this a reduced gravity model.
626  if (use_eos) then
627  !$OMP parallel do default(shared)
628  do j=jsq,jeq+1
629  if (use_p_atm) then
630  call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, &
631  tv%eqn_of_state, eosdom)
632  else
633  call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, &
634  tv%eqn_of_state, eosdom)
635  endif
636  do i=isq,ieq+1
637  dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * rho_in_situ(i)) * e(i,j,1)
638  enddo
639  enddo
640  else
641  !$OMP parallel do default(shared)
642  do j=jsq,jeq+1 ; do i=isq,ieq+1
643  dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * gv%Rlay(1)) * e(i,j,1)
644  enddo ; enddo
645  endif
646  endif
647  ! I have checked that rho_0 drops out and that the 1-layer case is right. RWH.
648 
649  ! If regridding is activated, do a linear reconstruction of salinity
650  ! and temperature across each layer. The subscripts 't' and 'b' refer
651  ! to top and bottom values within each layer (these are the only degrees
652  ! of freedeom needed to know the linear profile).
653  if ( use_ale ) then
654  if ( cs%Recon_Scheme == 1 ) then
655  call ts_plm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
656  elseif ( cs%Recon_Scheme == 2 ) then
657  call ts_ppm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
658  endif
659  endif
660 
661  ! Set the surface boundary conditions on pressure anomaly and its horizontal
662  ! integrals, assuming that the surface pressure anomaly varies linearly
663  ! in x and y.
664  if (use_p_atm) then
665  !$OMP parallel do default(shared)
666  do j=jsq,jeq+1 ; do i=isq,ieq+1
667  pa(i,j) = (rho_ref*gv%g_Earth)*e(i,j,1) + p_atm(i,j)
668  enddo ; enddo
669  else
670  !$OMP parallel do default(shared)
671  do j=jsq,jeq+1 ; do i=isq,ieq+1
672  pa(i,j) = (rho_ref*gv%g_Earth)*e(i,j,1)
673  enddo ; enddo
674  endif
675  !$OMP parallel do default(shared)
676  do j=js,je ; do i=isq,ieq
677  intx_pa(i,j) = 0.5*(pa(i,j) + pa(i+1,j))
678  enddo ; enddo
679  !$OMP parallel do default(shared)
680  do j=jsq,jeq ; do i=is,ie
681  inty_pa(i,j) = 0.5*(pa(i,j) + pa(i,j+1))
682  enddo ; enddo
683 
684  do k=1,nz
685  ! Calculate 4 integrals through the layer that are required in the
686  ! subsequent calculation.
687 
688  if (use_eos) then
689  ! The following routine computes the integrals that are needed to
690  ! calculate the pressure gradient force. Linear profiles for T and S are
691  ! assumed when regridding is activated. Otherwise, the previous version
692  ! is used, whereby densities within each layer are constant no matter
693  ! where the layers are located.
694  if ( use_ale ) then
695  if ( cs%Recon_Scheme == 1 ) then
696  call int_density_dz_generic_plm(k, tv, t_t, t_b, s_t, s_b, e, &
697  rho_ref, cs%Rho0, gv%g_Earth, dz_neglect, g%bathyT, &
698  g%HI, gv, tv%eqn_of_state, us, dpa, intz_dpa, intx_dpa, inty_dpa, &
699  usemasswghtinterp=cs%useMassWghtInterp)
700  elseif ( cs%Recon_Scheme == 2 ) then
701  call int_density_dz_generic_ppm(k, tv, t_t, t_b, s_t, s_b, e, &
702  rho_ref, cs%Rho0, gv%g_Earth, dz_neglect, g%bathyT, &
703  g%HI, gv, tv%eqn_of_state, us, dpa, intz_dpa, intx_dpa, inty_dpa, &
704  usemasswghtinterp=cs%useMassWghtInterp)
705  endif
706  else
707  call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,k), e(:,:,k+1), &
708  rho_ref, cs%Rho0, gv%g_Earth, g%HI, tv%eqn_of_state, us, dpa, &
709  intz_dpa, intx_dpa, inty_dpa, g%bathyT, dz_neglect, cs%useMassWghtInterp)
710  endif
711  !$OMP parallel do default(shared)
712  do j=jsq,jeq+1 ; do i=isq,ieq+1
713  intz_dpa(i,j) = intz_dpa(i,j)*gv%Z_to_H
714  enddo ; enddo
715  else
716  !$OMP parallel do default(shared)
717  do j=jsq,jeq+1 ; do i=isq,ieq+1
718  dz_geo(i,j) = gv%g_Earth * gv%H_to_Z*h(i,j,k)
719  dpa(i,j) = (gv%Rlay(k) - rho_ref) * dz_geo(i,j)
720  intz_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * dz_geo(i,j)*h(i,j,k)
721  enddo ; enddo
722  !$OMP parallel do default(shared)
723  do j=js,je ; do i=isq,ieq
724  intx_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i+1,j))
725  enddo ; enddo
726  !$OMP parallel do default(shared)
727  do j=jsq,jeq ; do i=is,ie
728  inty_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i,j+1))
729  enddo ; enddo
730  endif
731 
732  ! Compute pressure gradient in x direction
733  !$OMP parallel do default(shared)
734  do j=js,je ; do i=isq,ieq
735  pfu(i,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
736  (pa(i+1,j)*h(i+1,j,k) + intz_dpa(i+1,j))) + &
737  ((h(i+1,j,k) - h(i,j,k)) * intx_pa(i,j) - &
738  (e(i+1,j,k+1) - e(i,j,k+1)) * intx_dpa(i,j) * gv%Z_to_H)) * &
739  ((2.0*i_rho0*g%IdxCu(i,j)) / &
740  ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
741  intx_pa(i,j) = intx_pa(i,j) + intx_dpa(i,j)
742  enddo ; enddo
743  ! Compute pressure gradient in y direction
744  !$OMP parallel do default(shared)
745  do j=jsq,jeq ; do i=is,ie
746  pfv(i,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
747  (pa(i,j+1)*h(i,j+1,k) + intz_dpa(i,j+1))) + &
748  ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,j) - &
749  (e(i,j+1,k+1) - e(i,j,k+1)) * inty_dpa(i,j) * gv%Z_to_H)) * &
750  ((2.0*i_rho0*g%IdyCv(i,j)) / &
751  ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
752  inty_pa(i,j) = inty_pa(i,j) + inty_dpa(i,j)
753  enddo ; enddo
754  !$OMP parallel do default(shared)
755  do j=jsq,jeq+1 ; do i=isq,ieq+1
756  pa(i,j) = pa(i,j) + dpa(i,j)
757  enddo ; enddo
758  enddo
759 
760  if (cs%GFS_scale < 1.0) then
761  do k=1,nz
762  !$OMP parallel do default(shared)
763  do j=js,je ; do i=isq,ieq
764  pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
765  enddo ; enddo
766  !$OMP parallel do default(shared)
767  do j=jsq,jeq ; do i=is,ie
768  pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
769  enddo ; enddo
770  enddo
771  endif
772 
773  if (present(pbce)) then
774  call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce)
775  endif
776 
777  if (present(eta)) then
778  if (cs%tides) then
779  ! eta is the sea surface height relative to a time-invariant geoid, for comparison with
780  ! what is used for eta in btstep. See how e was calculated about 200 lines above.
781  !$OMP parallel do default(shared)
782  do j=jsq,jeq+1 ; do i=isq,ieq+1
783  eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
784  enddo ; enddo
785  else
786  !$OMP parallel do default(shared)
787  do j=jsq,jeq+1 ; do i=isq,ieq+1
788  eta(i,j) = e(i,j,1)*gv%Z_to_H
789  enddo ; enddo
790  endif
791  endif
792 
793  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
794  if (cs%id_tvar_sgs>0) call post_data(cs%id_tvar_sgs, tv%varT, cs%diag)
795 

◆ pressureforce_fv_end()

subroutine, public mom_pressureforce_fv::pressureforce_fv_end ( type(pressureforce_fv_cs), pointer  CS)

Deallocates the finite volume pressure gradient control structure.

Parameters
csFinite volume pressure control structure that will be deallocated in this subroutine.

Definition at line 880 of file MOM_PressureForce_FV.F90.

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

◆ pressureforce_fv_init()

subroutine, public mom_pressureforce_fv::pressureforce_fv_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_fv_cs), pointer  CS,
type(tidal_forcing_cs), optional, pointer  tides_CSp 
)

Initializes the finite volume pressure gradient control structure.

Parameters
[in]timeCurrent model time
[in]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]param_fileParameter file handles
[in,out]diagDiagnostics control structure
csFinite volume PGF control structure
tides_cspTides control structure

Definition at line 799 of file MOM_PressureForce_FV.F90.

800  type(time_type), target, intent(in) :: Time !< Current model time
801  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
802  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
803  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
804  type(param_file_type), intent(in) :: param_file !< Parameter file handles
805  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
806  type(PressureForce_FV_CS), pointer :: CS !< Finite volume PGF control structure
807  type(tidal_forcing_CS), optional, pointer :: tides_CSp !< Tides control structure
808  ! This include declares and sets the variable "version".
809 # include "version_variable.h"
810  character(len=40) :: mdl ! This module's name.
811  logical :: use_ALE
812 
813  if (associated(cs)) then
814  call mom_error(warning, "PressureForce_init called with an associated "// &
815  "control structure.")
816  return
817  else ; allocate(cs) ; endif
818 
819  cs%diag => diag ; cs%Time => time
820  if (present(tides_csp)) then
821  if (associated(tides_csp)) cs%tides_CSp => tides_csp
822  endif
823 
824  mdl = "MOM_PressureForce_FV"
825  call log_version(param_file, mdl, version, "")
826  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
827  "The mean ocean density used with BOUSSINESQ true to "//&
828  "calculate accelerations and the mass for conservation "//&
829  "properties, or with BOUSSINSEQ false to convert some "//&
830  "parameters from vertical units of m to kg m-2.", &
831  units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
832  call get_param(param_file, mdl, "TIDES", cs%tides, &
833  "If true, apply tidal momentum forcing.", default=.false.)
834  call get_param(param_file, "MOM", "USE_REGRIDDING", use_ale, &
835  "If True, use the ALE algorithm (regridding/remapping). "//&
836  "If False, use the layered isopycnal algorithm.", default=.false. )
837  call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", cs%useMassWghtInterp, &
838  "If true, use mass weighting when interpolating T/S for "//&
839  "integrals near the bathymetry in FV pressure gradient "//&
840  "calculations.", default=.false.)
841  call get_param(param_file, mdl, "RECONSTRUCT_FOR_PRESSURE", cs%reconstruct, &
842  "If True, use vertical reconstruction of T & S within "//&
843  "the integrals of the FV pressure gradient calculation. "//&
844  "If False, use the constant-by-layer algorithm. "//&
845  "The default is set by USE_REGRIDDING.", &
846  default=use_ale )
847  call get_param(param_file, mdl, "PRESSURE_RECONSTRUCTION_SCHEME", cs%Recon_Scheme, &
848  "Order of vertical reconstruction of T/S to use in the "//&
849  "integrals within the FV pressure gradient calculation.\n"//&
850  " 0: PCM or no reconstruction.\n"//&
851  " 1: PLM reconstruction.\n"//&
852  " 2: PPM reconstruction.", default=1)
853  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION_PRESSURE", cs%boundary_extrap, &
854  "If true, the reconstruction of T & S for pressure in "//&
855  "boundary cells is extrapolated, rather than using PCM "//&
856  "in these cells. If true, the same order polynomial is "//&
857  "used as is used for the interior cells.", default=.true.)
858  call get_param(param_file, mdl, "PGF_STANLEY_T2_DET_COEFF", cs%Stanley_T2_det_coeff, &
859  "The coefficient correlating SGS temperature variance with "// &
860  "the mean temperature gradient in the deterministic part of "// &
861  "the Stanley form of the Brankart correction. "// &
862  "Negative values disable the scheme.", units="nondim", default=-1.0)
863  if (cs%Stanley_T2_det_coeff>=0.) then
864  cs%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs_pgf', diag%axesTL, &
865  time, 'SGS temperature variance used in PGF', 'degC2')
866  endif
867  if (cs%tides) then
868  cs%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
869  time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=us%Z_to_m)
870  endif
871 
872  cs%GFS_scale = 1.0
873  if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
874 
875  call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale)
876 

◆ pressureforce_fv_nonbouss()

subroutine, public mom_pressureforce_fv::pressureforce_fv_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_fv_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]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thickness [H ~> kg/m2]
[in]tvThermodynamic variables
[out]pfuZonal acceleration [L T-2 ~> m s-2]
[out]pfvMeridional acceleration [L T-2 ~> m s-2]
csFinite volume PGF control structure
ale_cspALE control structure
p_atmThe pressure at the ice-ocean or atmosphere-ocean interface [R L2 T-2 ~> Pa].
[out]pbceThe baroclinic pressure anomaly in each layer due to eta anomalies [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
[out]etaThe bottom mass used to calculate PFu and PFv [H ~> m or kg m-2], with any tidal contributions or compressibility compensation.

Definition at line 76 of file MOM_PressureForce_FV.F90.

77  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
78  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
79  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
80  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> kg/m2]
81  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
82  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: PFu !< Zonal acceleration [L T-2 ~> m s-2]
83  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> m s-2]
84  type(PressureForce_FV_CS), pointer :: CS !< Finite volume PGF control structure
85  type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
86  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean
87  !! or atmosphere-ocean interface [R L2 T-2 ~> Pa].
88  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
89  !! anomaly in each layer due to eta anomalies
90  !! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
91  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
92  !! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
93  !! contributions or compressibility compensation.
94  ! Local variables
95  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure [R L2 T-2 ~> Pa].
96  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
97  T_tmp, & ! Temporary array of temperatures where layers that are lighter
98  ! than the mixed layer have the mixed layer's properties [degC].
99  s_tmp ! Temporary array of salinities where layers that are lighter
100  ! than the mixed layer have the mixed layer's properties [ppt].
101  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
102  S_t, & ! Top and bottom edge values for linear reconstructions
103  S_b, & ! of salinity within each layer [ppt].
104  T_t, & ! Top and bottom edge values for linear reconstructions
105  T_b ! of temperature within each layer [degC].
106  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
107  dza, & ! The change in geopotential anomaly between the top and bottom
108  ! of a layer [L2 T-2 ~> m2 s-2].
109  intp_dza ! The vertical integral in depth of the pressure anomaly less
110  ! the pressure anomaly at the top of the layer [R L4 Z-4 ~> Pa m2 s-2].
111  real, dimension(SZI_(G),SZJ_(G)) :: &
112  dp, & ! The (positive) change in pressure across a layer [R L2 Z-2 ~> Pa].
113  SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
114  e_tidal, & ! The bottom geopotential anomaly due to tidal forces from
115  ! astronomical sources and self-attraction and loading [Z ~> m].
116  dm, & ! The barotropic adjustment to the Montgomery potential to
117  ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
118  za ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the
119  ! interface atop a layer [L2 T-2 ~> m2 s-2].
120 
121  real, dimension(SZI_(G)) :: Rho_cv_BL ! The coordinate potential density in the deepest variable
122  ! density near-surface layer [R ~> kg m-3].
123  real, dimension(SZIB_(G),SZJ_(G)) :: &
124  intx_za ! The zonal integral of the geopotential anomaly along the
125  ! interface below a layer, divided by the grid spacing [L2 T-2 ~> m2 s-2].
126  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
127  intx_dza ! The change in intx_za through a layer [L2 T-2 ~> m2 s-2].
128  real, dimension(SZI_(G),SZJB_(G)) :: &
129  inty_za ! The meridional integral of the geopotential anomaly along the
130  ! interface below a layer, divided by the grid spacing [L2 T-2 ~> m2 s-2].
131  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
132  inty_dza ! The change in inty_za through a layer [L2 T-2 ~> m2 s-2].
133  real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
134  ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
135 
136  real :: dp_neglect ! A thickness that is so small it is usually lost
137  ! in roundoff and can be neglected [R L2 T-2 ~> Pa].
138  real :: I_gEarth ! The inverse of GV%g_Earth [L2 Z L-2 ~> s2 m-1]
139  real :: alpha_anom ! The in-situ specific volume, averaged over a
140  ! layer, less alpha_ref [R-1 ~> m3 kg-1].
141  logical :: use_p_atm ! If true, use the atmospheric pressure.
142  logical :: use_ALE ! If true, use an ALE pressure reconstruction.
143  logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
144  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
145 
146  real :: alpha_ref ! A reference specific volume [R-1 ~> m3 kg-1] that is used
147  ! to reduce the impact of truncation errors.
148  real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3].
149  real :: Pa_to_H ! A factor to convert from Pa to the thicknesss units (H) [H T2 R-1 L-2 ~> H Pa-1].
150  real :: H_to_RL2_T2 ! A factor to convert from thicknesss units (H) to pressure units [R L2 T-2 H-1 ~> Pa H-1].
151 ! real :: oneatm = 101325.0 ! 1 atm in [Pa] = [kg m-1 s-2]
152  real, parameter :: C1_6 = 1.0/6.0
153  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
154  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
155  integer :: i, j, k
156 
157  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
158  nkmb=gv%nk_rho_varies
159  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
160  eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)
161 
162  if (.not.associated(cs)) call mom_error(fatal, &
163  "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")
164  if (cs%Stanley_T2_det_coeff>=0.) call mom_error(fatal, &
165  "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//&
166  "implemented in non-Boussinesq mode.")
167 
168  use_p_atm = .false.
169  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
170  use_eos = associated(tv%eqn_of_state)
171  use_ale = .false.
172  if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
173 
174  h_to_rl2_t2 = gv%g_Earth*gv%H_to_RZ
175  dp_neglect = gv%g_Earth*gv%H_to_RZ * gv%H_subroundoff
176  alpha_ref = 1.0 / cs%Rho0
177  i_gearth = 1.0 / gv%g_Earth
178 
179  if (use_p_atm) then
180  !$OMP parallel do default(shared)
181  do j=jsq,jeq+1 ; do i=isq,ieq+1
182  p(i,j,1) = p_atm(i,j)
183  enddo ; enddo
184  else
185  !$OMP parallel do default(shared)
186  do j=jsq,jeq+1 ; do i=isq,ieq+1
187  p(i,j,1) = 0.0 ! or oneatm
188  enddo ; enddo
189  endif
190  !$OMP parallel do default(shared)
191  do j=jsq,jeq+1 ; do k=2,nz+1 ; do i=isq,ieq+1
192  p(i,j,k) = p(i,j,k-1) + h_to_rl2_t2 * h(i,j,k-1)
193  enddo ; enddo ; enddo
194 
195  if (use_eos) then
196  ! With a bulk mixed layer, replace the T & S of any layers that are
197  ! lighter than the the buffer layer with the properties of the buffer
198  ! layer. These layers will be massless anyway, and it avoids any
199  ! formal calculations with hydrostatically unstable profiles.
200  if (nkmb>0) then
201  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
202  tv_tmp%eqn_of_state => tv%eqn_of_state
203  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
204  !$OMP parallel do default(shared) private(Rho_cv_BL)
205  do j=jsq,jeq+1
206  do k=1,nkmb ; do i=isq,ieq+1
207  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
208  enddo ; enddo
209  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, rho_cv_bl(:), &
210  tv%eqn_of_state, eosdom)
211  do k=nkmb+1,nz ; do i=isq,ieq+1
212  if (gv%Rlay(k) < rho_cv_bl(i)) then
213  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
214  else
215  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
216  endif
217  enddo ; enddo
218  enddo
219  else
220  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
221  tv_tmp%eqn_of_state => tv%eqn_of_state
222  endif
223  endif
224 
225  ! If regridding is activated, do a linear reconstruction of salinity
226  ! and temperature across each layer. The subscripts 't' and 'b' refer
227  ! to top and bottom values within each layer (these are the only degrees
228  ! of freedeom needed to know the linear profile).
229  if ( use_ale ) then
230  if ( cs%Recon_Scheme == 1 ) then
231  call ts_plm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
232  elseif ( cs%Recon_Scheme == 2) then
233  call ts_ppm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
234  endif
235  endif
236 
237  !$OMP parallel do default(shared) private(alpha_anom,dp)
238  do k=1,nz
239  ! Calculate 4 integrals through the layer that are required in the
240  ! subsequent calculation.
241  if (use_eos) then
242  if ( use_ale ) then
243  if ( cs%Recon_Scheme == 1 ) then
244  call int_spec_vol_dp_generic_plm( t_t(:,:,k), t_b(:,:,k), s_t(:,:,k), s_b(:,:,k), &
245  p(:,:,k), p(:,:,k+1), alpha_ref, dp_neglect, p(:,:,nz+1), g%HI, &
246  tv%eqn_of_state, us, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), &
247  usemasswghtinterp=cs%useMassWghtInterp)
248  elseif ( cs%Recon_Scheme == 2 ) then
249  call mom_error(fatal, "PressureForce_FV_nonBouss: "//&
250  "int_spec_vol_dp_generic_ppm does not exist yet.")
251  ! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
252  ! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), &
253  ! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
254  ! intx_dza(:,:,k), inty_dza(:,:,k))
255  endif
256  else
257  call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,k), &
258  p(:,:,k+1), alpha_ref, g%HI, tv%eqn_of_state, &
259  us, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
260  inty_dza(:,:,k), bathyp=p(:,:,nz+1), dp_tiny=dp_neglect, &
261  usemasswghtinterp=cs%useMassWghtInterp)
262  endif
263  else
264  alpha_anom = 1.0 / gv%Rlay(k) - alpha_ref
265  do j=jsq,jeq+1 ; do i=isq,ieq+1
266  dp(i,j) = h_to_rl2_t2 * h(i,j,k)
267  dza(i,j,k) = alpha_anom * dp(i,j)
268  intp_dza(i,j,k) = 0.5 * alpha_anom * dp(i,j)**2
269  enddo ; enddo
270  do j=js,je ; do i=isq,ieq
271  intx_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i+1,j))
272  enddo ; enddo
273  do j=jsq,jeq ; do i=is,ie
274  inty_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i,j+1))
275  enddo ; enddo
276  endif
277  enddo
278 
279  ! The bottom geopotential anomaly is calculated first so that the increments
280  ! to the geopotential anomalies can be reused. Alternately, the surface
281  ! geopotential could be calculated directly with separate calls to
282  ! int_specific_vol_dp with alpha_ref=0, and the anomalies used going
283  ! downward, which would relieve the need for dza, intp_dza, intx_dza, and
284  ! inty_dza to be 3-D arrays.
285 
286  ! Sum vertically to determine the surface geopotential anomaly.
287  !$OMP parallel do default(shared)
288  do j=jsq,jeq+1
289  do i=isq,ieq+1
290  za(i,j) = alpha_ref*p(i,j,nz+1) - gv%g_Earth*g%bathyT(i,j)
291  enddo
292  do k=nz,1,-1 ; do i=isq,ieq+1
293  za(i,j) = za(i,j) + dza(i,j,k)
294  enddo ; enddo
295  enddo
296 
297  if (cs%tides) then
298  ! Find and add the tidal geopotential anomaly.
299  !$OMP parallel do default(shared)
300  do j=jsq,jeq+1 ; do i=isq,ieq+1
301  ssh(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * i_gearth
302  enddo ; enddo
303  call calc_tidal_forcing(cs%Time, ssh, e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
304  !$OMP parallel do default(shared)
305  do j=jsq,jeq+1 ; do i=isq,ieq+1
306  za(i,j) = za(i,j) - gv%g_Earth * e_tidal(i,j)
307  enddo ; enddo
308  endif
309 
310  if (cs%GFS_scale < 1.0) then
311  ! Adjust the Montgomery potential to make this a reduced gravity model.
312  if (use_eos) then
313  !$OMP parallel do default(shared) private(rho_in_situ)
314  do j=jsq,jeq+1
315  call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, &
316  tv%eqn_of_state, eosdom)
317 
318  do i=isq,ieq+1
319  dm(i,j) = (cs%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
320  enddo
321  enddo
322  else
323  !$OMP parallel do default(shared)
324  do j=jsq,jeq+1 ; do i=isq,ieq+1
325  dm(i,j) = (cs%GFS_scale - 1.0) * (p(i,j,1)*(1.0/gv%Rlay(1) - alpha_ref) + za(i,j))
326  enddo ; enddo
327  endif
328 ! else
329 ! do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; dM(i,j) = 0.0 ; enddo ; enddo
330  endif
331 
332  ! This order of integrating upward and then downward again is necessary with
333  ! a nonlinear equation of state, so that the surface geopotentials will go
334  ! linearly between the values at thickness points, but the bottom
335  ! geopotentials will not now be linear at the sub-grid-scale. Doing this
336  ! ensures no motion with flat isopycnals, even with a nonlinear equation of state.
337  !$OMP parallel do default(shared)
338  do j=js,je ; do i=isq,ieq
339  intx_za(i,j) = 0.5*(za(i,j) + za(i+1,j))
340  enddo ; enddo
341  !$OMP parallel do default(shared)
342  do j=jsq,jeq ; do i=is,ie
343  inty_za(i,j) = 0.5*(za(i,j) + za(i,j+1))
344  enddo ; enddo
345  do k=1,nz
346  ! These expressions for the acceleration have been carefully checked in
347  ! a set of idealized cases, and should be bug-free.
348  !$OMP parallel do default(shared)
349  do j=jsq,jeq+1 ; do i=isq,ieq+1
350  dp(i,j) = h_to_rl2_t2 * h(i,j,k)
351  za(i,j) = za(i,j) - dza(i,j,k)
352  enddo ; enddo
353  !$OMP parallel do default(shared)
354  do j=js,je ; do i=isq,ieq
355  intx_za(i,j) = intx_za(i,j) - intx_dza(i,j,k)
356  pfu(i,j,k) = ( ((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
357  (za(i+1,j)*dp(i+1,j) + intp_dza(i+1,j,k))) + &
358  ((dp(i+1,j) - dp(i,j)) * intx_za(i,j) - &
359  (p(i+1,j,k) - p(i,j,k)) * intx_dza(i,j,k)) ) * &
360  (2.0*g%IdxCu(i,j) / ((dp(i,j) + dp(i+1,j)) + dp_neglect))
361  enddo ; enddo
362  !$OMP parallel do default(shared)
363  do j=jsq,jeq ; do i=is,ie
364  inty_za(i,j) = inty_za(i,j) - inty_dza(i,j,k)
365  pfv(i,j,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
366  (za(i,j+1)*dp(i,j+1) + intp_dza(i,j+1,k))) + &
367  ((dp(i,j+1) - dp(i,j)) * inty_za(i,j) - &
368  (p(i,j+1,k) - p(i,j,k)) * inty_dza(i,j,k))) * &
369  (2.0*g%IdyCv(i,j) / ((dp(i,j) + dp(i,j+1)) + dp_neglect))
370  enddo ; enddo
371 
372  if (cs%GFS_scale < 1.0) then
373  ! Adjust the Montgomery potential to make this a reduced gravity model.
374  !$OMP parallel do default(shared)
375  do j=js,je ; do i=isq,ieq
376  pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
377  enddo ; enddo
378  !$OMP parallel do default(shared)
379  do j=jsq,jeq ; do i=is,ie
380  pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
381  enddo ; enddo
382  endif
383  enddo
384 
385  if (present(pbce)) then
386  call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce)
387  endif
388 
389  if (present(eta)) then
390  pa_to_h = 1.0 / (gv%g_Earth * gv%H_to_RZ)
391  if (use_p_atm) then
392  !$OMP parallel do default(shared)
393  do j=jsq,jeq+1 ; do i=isq,ieq+1
394  eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h ! eta has the same units as h.
395  enddo ; enddo
396  else
397  !$OMP parallel do default(shared)
398  do j=jsq,jeq+1 ; do i=isq,ieq+1
399  eta(i,j) = p(i,j,nz+1)*pa_to_h ! eta has the same units as h.
400  enddo ; enddo
401  endif
402  endif
403 
404  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
405