MOM6
Data Types | Functions/Subroutines
mom_pressureforce_blk_afv Module Reference

Detailed Description

Analytically integrated finite volume pressure gradient.

Data Types

type  pressureforce_blk_afv_cs
 Finite volume pressure gradient control structure. More...
 

Functions/Subroutines

subroutine, public pressureforce_blk_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_blk_afv_nonbouss (h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
 Non-Boussinesq analytically-integrated finite volume form of pressure gradient. More...
 
subroutine, public pressureforce_blk_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_blk_afv_init (Time, G, GV, US, param_file, diag, CS, tides_CSp)
 Initializes the finite volume pressure gradient control structure. More...
 
subroutine, public pressureforce_blk_afv_end (CS)
 Deallocates the finite volume pressure gradient control structure. More...
 

Function/Subroutine Documentation

◆ pressureforce_blk_afv()

subroutine, public mom_pressureforce_blk_afv::pressureforce_blk_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_blk_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]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thickness [H ~> m or kg m-2]
[in,out]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 [Pa].
[out]pbceThe 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]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 67 of file MOM_PressureForce_blocked_AFV.F90.

67  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
68  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
69  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
70  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
71  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variables
72  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: PFu !< Zonal acceleration [L T-2 ~> m s-2]
73  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> m s-2]
74  type(PressureForce_blk_AFV_CS), pointer :: CS !< Finite volume PGF control structure
75  type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
76  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean
77  !! or atmosphere-ocean interface [Pa].
78  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
79  !! anomaly in each layer due to eta anomalies
80  !! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
81  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
82  !! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
83  !! contributions or compressibility compensation.
84 
85  if (gv%Boussinesq) then
86  call pressureforce_blk_afv_bouss(h, tv, pfu, pfv, g, gv, us, cs, ale_csp, p_atm, pbce, eta)
87  else
88  call pressureforce_blk_afv_nonbouss(h, tv, pfu, pfv, g, gv, us, cs, p_atm, pbce, eta)
89  endif
90 

References pressureforce_blk_afv_bouss(), and pressureforce_blk_afv_nonbouss().

Here is the call graph for this function:

◆ pressureforce_blk_afv_bouss()

subroutine, public mom_pressureforce_blk_afv::pressureforce_blk_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_blk_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, making the Boussinesq approximation. This version uses code-blocking for threads.

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 or kg m-2]
[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 [Pa].
[out]pbceThe 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]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 427 of file MOM_PressureForce_blocked_AFV.F90.

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

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

Here is the call graph for this function:
Here is the caller graph for this function:

◆ pressureforce_blk_afv_end()

subroutine, public mom_pressureforce_blk_afv::pressureforce_blk_afv_end ( type(pressureforce_blk_afv_cs), pointer  CS)

Deallocates the finite volume pressure gradient control structure.

Parameters
csBlocked AFV pressure control structure that will be deallocated in this subroutine.

Definition at line 856 of file MOM_PressureForce_blocked_AFV.F90.

856  type(PressureForce_blk_AFV_CS), pointer :: CS !< Blocked AFV pressure control structure that
857  !! will be deallocated in this subroutine.
858  if (associated(cs)) deallocate(cs)

◆ pressureforce_blk_afv_init()

subroutine, public mom_pressureforce_blk_afv::pressureforce_blk_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_blk_afv_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 783 of file MOM_PressureForce_blocked_AFV.F90.

783  type(time_type), target, intent(in) :: Time !< Current model time
784  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
785  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
786  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
787  type(param_file_type), intent(in) :: param_file !< Parameter file handles
788  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
789  type(PressureForce_blk_AFV_CS), pointer :: CS !< Finite volume PGF control structure
790  type(tidal_forcing_CS), optional, pointer :: tides_CSp !< Tides control structure
791 ! This include declares and sets the variable "version".
792 #include "version_variable.h"
793  character(len=40) :: mdl ! This module's name.
794  logical :: use_ALE
795 
796  if (associated(cs)) then
797  call mom_error(warning, "PressureForce_init called with an associated "// &
798  "control structure.")
799  return
800  else ; allocate(cs) ; endif
801 
802  cs%diag => diag ; cs%Time => time
803  if (present(tides_csp)) then
804  if (associated(tides_csp)) cs%tides_CSp => tides_csp
805  endif
806 
807  mdl = "MOM_PressureForce_blk_AFV"
808  call log_version(param_file, mdl, version, "")
809  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
810  "The mean ocean density used with BOUSSINESQ true to "//&
811  "calculate accelerations and the mass for conservation "//&
812  "properties, or with BOUSSINSEQ false to convert some "//&
813  "parameters from vertical units of m to kg m-2.", &
814  units="kg m-3", default=1035.0)
815  call get_param(param_file, mdl, "TIDES", cs%tides, &
816  "If true, apply tidal momentum forcing.", default=.false.)
817  call get_param(param_file, "MOM", "USE_REGRIDDING", use_ale, &
818  "If True, use the ALE algorithm (regridding/remapping). "//&
819  "If False, use the layered isopycnal algorithm.", default=.false. )
820  call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", cs%useMassWghtInterp, &
821  "If true, use mass weighting when interpolating T/S for "//&
822  "integrals near the bathymetry in AFV pressure gradient "//&
823  "calculations.", default=.false.)
824  call get_param(param_file, mdl, "RECONSTRUCT_FOR_PRESSURE", cs%reconstruct, &
825  "If True, use vertical reconstruction of T & S within "//&
826  "the integrals of the FV pressure gradient calculation. "//&
827  "If False, use the constant-by-layer algorithm. "//&
828  "The default is set by USE_REGRIDDING.", &
829  default=use_ale )
830  call get_param(param_file, mdl, "PRESSURE_RECONSTRUCTION_SCHEME", cs%Recon_Scheme, &
831  "Order of vertical reconstruction of T/S to use in the "//&
832  "integrals within the FV pressure gradient calculation.\n"//&
833  " 0: PCM or no reconstruction.\n"//&
834  " 1: PLM reconstruction.\n"//&
835  " 2: PPM reconstruction.", default=1)
836  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION_PRESSURE", cs%boundary_extrap, &
837  "If true, the reconstruction of T & S for pressure in "//&
838  "boundary cells is extrapolated, rather than using PCM "//&
839  "in these cells. If true, the same order polynomial is "//&
840  "used as is used for the interior cells.", default=.true.)
841 
842  if (cs%tides) then
843  cs%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
844  time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=us%Z_to_m)
845  endif
846 
847  cs%GFS_scale = 1.0
848  if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
849 
850  call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale)
851 

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ pressureforce_blk_afv_nonbouss()

subroutine, public mom_pressureforce_blk_afv::pressureforce_blk_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_blk_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. This version uses code-blocking for threads.

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 m-2]
[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
p_atmThe pressure at the ice-ocean or atmosphere-ocean interface [Pa].
[out]pbceThe 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]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 103 of file MOM_PressureForce_blocked_AFV.F90.

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

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

Referenced by pressureforce_blk_afv().

Here is the call graph for this function:
Here is the caller graph for this function: