MOM6
MOM_PressureForce_Montgomery.F90
Go to the documentation of this file.
1 !> Provides the Montgomery potential form of pressure gradient
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
7 use mom_diag_mediator, only : safe_alloc_ptr, diag_ctrl, time_type
8 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
10 use mom_grid, only : ocean_grid_type
17 
18 implicit none ; private
19 
20 #include <MOM_memory.h>
21 
24 
25 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
26 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
27 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
28 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
29 
30 !> Control structure for the Montgomery potential form of pressure gradient
31 type, public :: pressureforce_mont_cs ; private
32  logical :: tides !< If true, apply tidal momentum forcing.
33  real :: rho0 !< The density used in the Boussinesq
34  !! approximation [kg m-3].
35  real :: rho_atm !< The assumed atmospheric density [kg m-3].
36  !! By default, Rho_atm is 0.
37  real :: gfs_scale !< Ratio between gravity applied to top interface and the
38  !! gravitational acceleration of the planet [nondim].
39  !! Usually this ratio is 1.
40  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
41  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate
42  !! the timing of diagnostic output.
43  real, pointer :: pfu_bc(:,:,:) => null() !< Zonal accelerations due to pressure gradients
44  !! deriving from density gradients within layers [L T-2 ~> m s-2].
45  real, pointer :: pfv_bc(:,:,:) => null() !< Meridional accelerations due to pressure gradients
46  !! deriving from density gradients within layers [L T-2 ~> m s-2].
47  !>@{ Diagnostic IDs
48  integer :: id_pfu_bc = -1, id_pfv_bc = -1, id_e_tidal = -1
49  !!@}
50  type(tidal_forcing_cs), pointer :: tides_csp => null() !< The tidal forcing control structure
51 end type pressureforce_mont_cs
52 
53 contains
54 
55 !> \brief Non-Boussinesq Montgomery-potential form of pressure gradient
56 !!
57 !! Determines the acceleration due to pressure forces in a
58 !! non-Boussinesq fluid using the compressibility compensated (if appropriate)
59 !! Montgomery-potential form described in Hallberg (Ocean Mod., 2005).
60 !!
61 !! To work, the following fields must be set outside of the usual (is:ie,js:je)
62 !! range before this subroutine is called:
63 !! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
64 subroutine pressureforce_mont_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
65  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure.
66  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
67  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
68  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, [H ~> kg m-2].
69  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables.
70  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: pfu !< Zonal acceleration due to pressure gradients
71  !! (equal to -dM/dx) [L T-2 ~> m s-2].
72  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration due to pressure gradients
73  !! (equal to -dM/dy) [L T-2 ~> m s-2].
74  type(pressureforce_mont_cs), pointer :: cs !< Control structure for Montgomery potential PGF
75  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean or
76  !! atmosphere-ocean [Pa].
77  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
78  optional, intent(out) :: pbce !< The baroclinic pressure anomaly in
79  !! each layer due to free surface height anomalies,
80  !! [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2].
81  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< Free surface height [H ~> kg m-1].
82 
83  ! Local variables
84  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
85  m, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2].
86  alpha_star, & ! Compression adjusted specific volume [R-1 ~> m3 kg-1].
87  dz_geo ! The change in geopotential across a layer [m2 s-2].
88  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure [Pa].
89  ! p may be adjusted (with a nonlinear equation of state) so that
90  ! its derivative compensates for the adiabatic compressibility
91  ! in seawater, but p will still be close to the pressure.
92  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
93  t_tmp, & ! Temporary array of temperatures where layers that are lighter
94  ! than the mixed layer have the mixed layer's properties [degC].
95  s_tmp ! Temporary array of salinities where layers that are lighter
96  ! than the mixed layer have the mixed layer's properties [ppt].
97 
98  real, dimension(SZI_(G)) :: rho_cv_bl ! The coordinate potential density in the
99  ! deepest variable density near-surface layer [R ~> kg m-3].
100 
101  real, dimension(SZI_(G),SZJ_(G)) :: &
102  dm, & ! A barotropic correction to the Montgomery potentials to
103  ! enable the use of a reduced gravity form of the equations
104  ! [m2 s-2].
105  dp_star, & ! Layer thickness after compensation for compressibility [Pa].
106  ssh, & ! The sea surface height anomaly, in depth units [Z ~> m].
107  e_tidal, & ! Bottom geopotential anomaly due to tidal forces from
108  ! astronomical sources and self-attraction and loading [Z ~> m].
109  geopot_bot ! Bottom geopotential relative to time-mean sea level,
110  ! including any tidal contributions [L2 T-2 ~> m2 s-2].
111  real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate
112  ! density [Pa] (usually 2e7 Pa = 2000 dbar).
113  real :: rho_in_situ(szi_(g)) !In-situ density of a layer [R ~> kg m-3].
114  real :: pfu_bc, pfv_bc ! The pressure gradient force due to along-layer
115  ! compensated density gradients [L T-2 ~> m s-2]
116  real :: dp_neglect ! A thickness that is so small it is usually lost
117  ! in roundoff and can be neglected [Pa].
118  logical :: use_p_atm ! If true, use the atmospheric pressure.
119  logical :: use_eos ! If true, density is calculated from T & S using
120  ! an equation of state.
121  logical :: is_split ! A flag indicating whether the pressure
122  ! gradient terms are to be split into
123  ! barotropic and baroclinic pieces.
124  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
125 
126  real :: i_gearth ! The inverse of g_Earth [s2 Z m-2 ~> s2 m-1]
127 ! real :: dalpha
128  real :: pa_to_p_dyn ! A conversion factor from Pa (= kg m-1 s-2) to the units of
129  ! dynamic pressure (R L2 T-2) [ R L2 T-2 m s2 kg-1 ~> nondim]
130  real :: pa_to_h ! A factor to convert from Pa to the thicknesss units (H).
131  real :: alpha_lay(szk_(g)) ! The specific volume of each layer [R-1 ~> m3 kg-1].
132  real :: dalpha_int(szk_(g)+1) ! The change in specific volume across each
133  ! interface [R-1 ~> m3 kg-1].
134  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
135  integer :: i, j, k
136  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
137  nkmb=gv%nk_rho_varies
138  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
139 
140  use_p_atm = .false.
141  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
142  is_split = .false. ; if (present(pbce)) is_split = .true.
143  use_eos = associated(tv%eqn_of_state)
144 
145  if (.not.associated(cs)) call mom_error(fatal, &
146  "MOM_PressureForce_Mont: Module must be initialized before it is used.")
147  if (use_eos) then
148  if (query_compressible(tv%eqn_of_state)) call mom_error(fatal, &
149  "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//&
150  "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
151  endif
152 
153  pa_to_p_dyn = us%kg_m3_to_R * us%m_s_to_L_T**2
154  i_gearth = 1.0 / (us%L_T_to_m_s**2 * gv%g_Earth)
155  dp_neglect = gv%H_to_Pa * gv%H_subroundoff
156  do k=1,nz ; alpha_lay(k) = 1.0 / (gv%Rlay(k)) ; enddo
157  do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ; enddo
158 
159  if (use_p_atm) then
160  !$OMP parallel do default(shared)
161  do j=jsq,jeq+1 ; do i=isq,ieq+1 ; p(i,j,1) = p_atm(i,j) ; enddo ; enddo
162  else
163  !$OMP parallel do default(shared)
164  do j=jsq,jeq+1 ; do i=isq,ieq+1 ; p(i,j,1) = 0.0 ; enddo ; enddo
165  endif
166  !$OMP parallel do default(shared)
167  do j=jsq,jeq+1 ; do k=1,nz ; do i=isq,ieq+1
168  p(i,j,k+1) = p(i,j,k) + gv%H_to_Pa * h(i,j,k)
169  enddo ; enddo ; enddo
170 
171  if (present(eta)) then
172  pa_to_h = 1.0 / gv%H_to_Pa
173  if (use_p_atm) then
174  !$OMP parallel do default(shared)
175  do j=jsq,jeq+1 ; do i=isq,ieq+1
176  eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h ! eta has the same units as h.
177  enddo ; enddo
178  else
179  !$OMP parallel do default(shared)
180  do j=jsq,jeq+1 ; do i=isq,ieq+1
181  eta(i,j) = p(i,j,nz+1)*pa_to_h ! eta has the same units as h.
182  enddo ; enddo
183  endif
184  endif
185 
186  if (cs%tides) then
187  ! Determine the sea surface height anomalies, to enable the calculation
188  ! of self-attraction and loading.
189  !$OMP parallel do default(shared)
190  do j=jsq,jeq+1 ; do i=isq,ieq+1
191  ssh(i,j) = -g%bathyT(i,j)
192  enddo ; enddo
193  if (use_eos) then
194  !$OMP parallel do default(shared)
195  do k=1,nz
196  call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
197  0.0, g%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1)
198  enddo
199  !$OMP parallel do default(shared)
200  do j=jsq,jeq+1 ; do k=1,nz; do i=isq,ieq+1
201  ssh(i,j) = ssh(i,j) + i_gearth * dz_geo(i,j,k)
202  enddo ; enddo ; enddo
203  else
204  !$OMP parallel do default(shared)
205  do j=jsq,jeq+1 ; do k=1,nz ; do i=isq,ieq+1
206  ssh(i,j) = ssh(i,j) + gv%H_to_RZ * h(i,j,k) * alpha_lay(k)
207  enddo ; enddo ; enddo
208  endif
209 
210  call calc_tidal_forcing(cs%Time, ssh, e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
211  !$OMP parallel do default(shared)
212  do j=jsq,jeq+1 ; do i=isq,ieq+1
213  geopot_bot(i,j) = -gv%g_Earth*(e_tidal(i,j) + g%bathyT(i,j))
214  enddo ; enddo
215  else
216  !$OMP parallel do default(shared)
217  do j=jsq,jeq+1 ; do i=isq,ieq+1
218  geopot_bot(i,j) = -gv%g_Earth*g%bathyT(i,j)
219  enddo ; enddo
220  endif
221 
222  if (use_eos) then
223  ! Calculate in-situ specific volumes (alpha_star).
224 
225  ! With a bulk mixed layer, replace the T & S of any layers that are
226  ! lighter than the the buffer layer with the properties of the buffer
227  ! layer. These layers will be massless anyway, and it avoids any
228  ! formal calculations with hydrostatically unstable profiles.
229  if (nkmb>0) then
230  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
231  tv_tmp%eqn_of_state => tv%eqn_of_state
232  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
233  !$OMP parallel do default(shared) private(Rho_cv_BL)
234  do j=jsq,jeq+1
235  do k=1,nkmb ; do i=isq,ieq+1
236  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
237  enddo ; enddo
238  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
239  rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
240  do k=nkmb+1,nz ; do i=isq,ieq+1
241  if (gv%Rlay(k) < rho_cv_bl(i)) then
242  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
243  else
244  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
245  endif
246  enddo ; enddo
247  enddo
248  else
249  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
250  tv_tmp%eqn_of_state => tv%eqn_of_state
251  do i=isq,ieq+1 ; p_ref(i) = 0 ; enddo
252  endif
253  !$OMP parallel do default(shared) private(rho_in_situ)
254  do k=1,nz ; do j=jsq,jeq+1
255  call calculate_density(tv_tmp%T(:,j,k),tv_tmp%S(:,j,k),p_ref, &
256  rho_in_situ,isq,ieq-isq+2,tv%eqn_of_state, scale=us%kg_m3_to_R)
257  do i=isq,ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo
258  enddo ; enddo
259  endif ! use_EOS
260 
261  if (use_eos) then
262  !$OMP parallel do default(shared)
263  do j=jsq,jeq+1
264  do i=isq,ieq+1
265  m(i,j,nz) = geopot_bot(i,j) + pa_to_p_dyn*p(i,j,nz+1) * alpha_star(i,j,nz)
266  enddo
267  do k=nz-1,1,-1 ; do i=isq,ieq+1
268  m(i,j,k) = m(i,j,k+1) + pa_to_p_dyn*p(i,j,k+1) * (alpha_star(i,j,k) - alpha_star(i,j,k+1))
269  enddo ; enddo
270  enddo
271  else ! not use_EOS
272  !$OMP parallel do default(shared)
273  do j=jsq,jeq+1
274  do i=isq,ieq+1
275  m(i,j,nz) = geopot_bot(i,j) + pa_to_p_dyn*p(i,j,nz+1) * alpha_lay(nz)
276  enddo
277  do k=nz-1,1,-1 ; do i=isq,ieq+1
278  m(i,j,k) = m(i,j,k+1) + pa_to_p_dyn*p(i,j,k+1) * dalpha_int(k+1)
279  enddo ; enddo
280  enddo
281  endif ! use_EOS
282 
283  if (cs%GFS_scale < 1.0) then
284  ! Adjust the Montgomery potential to make this a reduced gravity model.
285  !$OMP parallel do default(shared)
286  do j=jsq,jeq+1 ; do i=isq,ieq+1
287  dm(i,j) = (cs%GFS_scale - 1.0) * m(i,j,1)
288  enddo ; enddo
289  !$OMP parallel do default(shared)
290  do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1
291  m(i,j,k) = m(i,j,k) + dm(i,j)
292  enddo ; enddo ; enddo
293 
294  ! Could instead do the following, to avoid taking small differences
295  ! of large numbers...
296 ! do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
297 ! M(i,j,1) = CS%GFS_scale * M(i,j,1)
298 ! enddo ; enddo
299 ! if (use_EOS) then
300 ! do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
301 ! M(i,j,k) = M(i,j,k-1) - Pa_to_p_dyn*p(i,j,K) * (alpha_star(i,j,k-1) - alpha_star(i,j,k))
302 ! enddo ; enddo ; enddo
303 ! else ! not use_EOS
304 ! do k=2,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
305 ! M(i,j,k) = M(i,j,k-1) - Pa_to_p_dyn*p(i,j,K) * dalpha_int(K)
306 ! enddo ; enddo ; enddo
307 ! endif ! use_EOS
308 
309  endif
310 
311  ! Note that ddM/dPb = alpha_star(i,j,1)
312  if (present(pbce)) then
313  call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce, alpha_star)
314  endif
315 
316 ! Calculate the pressure force. On a Cartesian grid,
317 ! PFu = - dM/dx and PFv = - dM/dy.
318  if (use_eos) then
319  !$OMP parallel do default(shared) private(dp_star,PFu_bc,PFv_bc)
320  do k=1,nz
321  do j=jsq,jeq+1 ; do i=isq,ieq+1
322  dp_star(i,j) = (p(i,j,k+1) - p(i,j,k)) + dp_neglect
323  enddo ; enddo
324  do j=js,je ; do i=isq,ieq
325  ! PFu_bc = p* grad alpha*
326  pfu_bc = (alpha_star(i+1,j,k) - alpha_star(i,j,k)) * (g%IdxCu(i,j) * pa_to_p_dyn * &
327  ((dp_star(i,j) * dp_star(i+1,j) + (p(i,j,k) * dp_star(i+1,j) + &
328  p(i+1,j,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i+1,j))))
329  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
330  if (associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
331  enddo ; enddo
332  do j=jsq,jeq ; do i=is,ie
333  pfv_bc = (alpha_star(i,j+1,k) - alpha_star(i,j,k)) * (g%IdyCv(i,j) * pa_to_p_dyn * &
334  ((dp_star(i,j) * dp_star(i,j+1) + (p(i,j,k) * dp_star(i,j+1) + &
335  p(i,j+1,k) * dp_star(i,j))) / (dp_star(i,j) + dp_star(i,j+1))))
336  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
337  if (associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
338  enddo ; enddo
339  enddo ! k-loop
340  else ! .not. use_EOS
341  !$OMP parallel do default(shared)
342  do k=1,nz
343  do j=js,je ; do i=isq,ieq
344  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
345  enddo ; enddo
346  do j=jsq,jeq ; do i=is,ie
347  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
348  enddo ; enddo
349  enddo
350  endif ! use_EOS
351 
352  if (cs%id_PFu_bc>0) call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
353  if (cs%id_PFv_bc>0) call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
354  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
355 
356 end subroutine pressureforce_mont_nonbouss
357 
358 !> \brief Boussinesq Montgomery-potential form of pressure gradient
359 !!
360 !! Determines the acceleration due to pressure forces.
361 !!
362 !! To work, the following fields must be set outside of the usual (is:ie,js:je)
363 !! range before this subroutine is called:
364 !! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
365 subroutine pressureforce_mont_bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
366  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure.
367  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
368  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
369  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m].
370  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables.
371  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: pfu !< Zonal acceleration due to pressure gradients
372  !! (equal to -dM/dx) [L T-2 ~> m s-2].
373  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration due to pressure gradients
374  !! (equal to -dM/dy) [L T-2 ~> m s2].
375  type(pressureforce_mont_cs), pointer :: cs !< Control structure for Montgomery potential PGF
376  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean or
377  !! atmosphere-ocean [Pa].
378  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure anomaly in
379  !! each layer due to free surface height anomalies
380  !! [L2 T-2 H-1 ~> m s-2].
381  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< Free surface height [H ~> m].
382  ! Local variables
383  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
384  m, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2].
385  rho_star ! In-situ density divided by the derivative with depth of the
386  ! corrected e times (G_Earth/Rho0) [m2 Z-1 s-2 ~> m s-2].
387  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in m.
388  ! e may be adjusted (with a nonlinear equation of state) so that
389  ! its derivative compensates for the adiabatic compressibility
390  ! in seawater, but e will still be close to the interface depth.
391  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
392  t_tmp, & ! Temporary array of temperatures where layers that are lighter
393  ! than the mixed layer have the mixed layer's properties [degC].
394  s_tmp ! Temporary array of salinities where layers that are lighter
395  ! than the mixed layer have the mixed layer's properties [ppt].
396 
397  real :: rho_cv_bl(szi_(g)) ! The coordinate potential density in
398  ! the deepest variable density near-surface layer [R ~> kg m-3].
399  real :: h_star(szi_(g),szj_(g)) ! Layer thickness after compensation
400  ! for compressibility [Z ~> m].
401  real :: e_tidal(szi_(g),szj_(g)) ! Bottom geopotential anomaly due to tidal
402  ! forces from astronomical sources and self-
403  ! attraction and loading, in depth units [Z ~> m].
404  real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate
405  ! density [Pa] (usually 2e7 Pa = 2000 dbar).
406  real :: i_rho0 ! 1/Rho0 [m3 kg-1].
407  real :: g_rho0 ! G_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
408  real :: pfu_bc, pfv_bc ! The pressure gradient force due to along-layer
409  ! compensated density gradients [L T-2 ~> m s-2]
410  real :: h_neglect ! A thickness that is so small it is usually lost
411  ! in roundoff and can be neglected [Z ~> m].
412  logical :: use_p_atm ! If true, use the atmospheric pressure.
413  logical :: use_eos ! If true, density is calculated from T & S using
414  ! an equation of state.
415  logical :: is_split ! A flag indicating whether the pressure
416  ! gradient terms are to be split into
417  ! barotropic and baroclinic pieces.
418  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
419  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
420  integer :: i, j, k
421 
422  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
423  nkmb=gv%nk_rho_varies
424  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
425 
426  use_p_atm = .false.
427  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
428  is_split = .false. ; if (present(pbce)) is_split = .true.
429  use_eos = associated(tv%eqn_of_state)
430 
431  if (.not.associated(cs)) call mom_error(fatal, &
432  "MOM_PressureForce_Mont: Module must be initialized before it is used.")
433  if (use_eos) then
434  if (query_compressible(tv%eqn_of_state)) call mom_error(fatal, &
435  "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//&
436  "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.")
437  endif
438 
439  h_neglect = gv%H_subroundoff * gv%H_to_Z
440  i_rho0 = 1.0/cs%Rho0
441  g_rho0 = gv%g_Earth / gv%Rho0
442 
443  if (cs%tides) then
444  ! Determine the surface height anomaly for calculating self attraction
445  ! and loading. This should really be based on bottom pressure anomalies,
446  ! but that is not yet implemented, and the current form is correct for
447  ! barotropic tides.
448  !$OMP parallel do default(shared)
449  do j=jsq,jeq+1
450  do i=isq,ieq+1 ; e(i,j,1) = -g%bathyT(i,j) ; enddo
451  do k=1,nz ; do i=isq,ieq+1
452  e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
453  enddo ; enddo
454  enddo
455  call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
456  endif
457 
458 ! Here layer interface heights, e, are calculated.
459  if (cs%tides) then
460  !$OMP parallel do default(shared)
461  do j=jsq,jeq+1 ; do i=isq,ieq+1
462  e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
463  enddo ; enddo
464  else
465  !$OMP parallel do default(shared)
466  do j=jsq,jeq+1 ; do i=isq,ieq+1
467  e(i,j,nz+1) = -g%bathyT(i,j)
468  enddo ; enddo
469  endif
470  !$OMP parallel do default(shared)
471  do j=jsq,jeq+1 ; do k=nz,1,-1 ; do i=isq,ieq+1
472  e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
473  enddo ; enddo ; enddo
474 
475  if (use_eos) then
476 ! Calculate in-situ densities (rho_star).
477 
478 ! With a bulk mixed layer, replace the T & S of any layers that are
479 ! lighter than the the buffer layer with the properties of the buffer
480 ! layer. These layers will be massless anyway, and it avoids any
481 ! formal calculations with hydrostatically unstable profiles.
482 
483  if (nkmb>0) then
484  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
485  tv_tmp%eqn_of_state => tv%eqn_of_state
486 
487  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
488  !$OMP parallel do default(shared) private(Rho_cv_BL)
489  do j=jsq,jeq+1
490  do k=1,nkmb ; do i=isq,ieq+1
491  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
492  enddo ; enddo
493  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
494  rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
495 
496  do k=nkmb+1,nz ; do i=isq,ieq+1
497  if (gv%Rlay(k) < rho_cv_bl(i)) then
498  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
499  else
500  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
501  endif
502  enddo ; enddo
503  enddo
504  else
505  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
506  tv_tmp%eqn_of_state => tv%eqn_of_state
507  do i=isq,ieq+1 ; p_ref(i) = 0.0 ; enddo
508  endif
509 
510  ! This no longer includes any pressure dependency, since this routine
511  ! will come down with a fatal error if there is any compressibility.
512  !$OMP parallel do default(shared)
513  do k=1,nz+1 ; do j=jsq,jeq+1
514  call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
515  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R*g_rho0)
516  enddo ; enddo
517  endif ! use_EOS
518 
519 ! Here the layer Montgomery potentials, M, are calculated.
520  if (use_eos) then
521  !$OMP parallel do default(shared)
522  do j=jsq,jeq+1
523  do i=isq,ieq+1
524  m(i,j,1) = cs%GFS_scale * (rho_star(i,j,1) * e(i,j,1))
525  if (use_p_atm) m(i,j,1) = m(i,j,1) + us%m_s_to_L_T**2*p_atm(i,j) * i_rho0
526  enddo
527  do k=2,nz ; do i=isq,ieq+1
528  m(i,j,k) = m(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,k)
529  enddo ; enddo
530  enddo
531  else ! not use_EOS
532  !$OMP parallel do default(shared)
533  do j=jsq,jeq+1
534  do i=isq,ieq+1
535  m(i,j,1) = gv%g_prime(1) * e(i,j,1)
536  if (use_p_atm) m(i,j,1) = m(i,j,1) + us%m_s_to_L_T**2*p_atm(i,j) * i_rho0
537  enddo
538  do k=2,nz ; do i=isq,ieq+1
539  m(i,j,k) = m(i,j,k-1) + gv%g_prime(k) * e(i,j,k)
540  enddo ; enddo
541  enddo
542  endif ! use_EOS
543 
544  if (present(pbce)) then
545  call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce, rho_star)
546  endif
547 
548 ! Calculate the pressure force. On a Cartesian grid,
549 ! PFu = - dM/dx and PFv = - dM/dy.
550  if (use_eos) then
551  !$OMP parallel do default(shared) private(h_star,PFu_bc,PFv_bc)
552  do k=1,nz
553  do j=jsq,jeq+1 ; do i=isq,ieq+1
554  h_star(i,j) = (e(i,j,k) - e(i,j,k+1)) + h_neglect
555  enddo ; enddo
556  do j=js,je ; do i=isq,ieq
557  pfu_bc = -1.0*(rho_star(i+1,j,k) - rho_star(i,j,k)) * (g%IdxCu(i,j) * &
558  ((h_star(i,j) * h_star(i+1,j) - (e(i,j,k) * h_star(i+1,j) + &
559  e(i+1,j,k) * h_star(i,j))) / (h_star(i,j) + h_star(i+1,j))))
560  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j) + pfu_bc
561  if (associated(cs%PFu_bc)) cs%PFu_bc(i,j,k) = pfu_bc
562  enddo ; enddo
563  do j=jsq,jeq ; do i=is,ie
564  pfv_bc = -1.0*(rho_star(i,j+1,k) - rho_star(i,j,k)) * (g%IdyCv(i,j) * &
565  ((h_star(i,j) * h_star(i,j+1) - (e(i,j,k) * h_star(i,j+1) + &
566  e(i,j+1,k) * h_star(i,j))) / (h_star(i,j) + h_star(i,j+1))))
567  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j) + pfv_bc
568  if (associated(cs%PFv_bc)) cs%PFv_bc(i,j,k) = pfv_bc
569  enddo ; enddo
570  enddo ! k-loop
571  else ! .not. use_EOS
572  !$OMP parallel do default(shared)
573  do k=1,nz
574  do j=js,je ; do i=isq,ieq
575  pfu(i,j,k) = -(m(i+1,j,k) - m(i,j,k)) * g%IdxCu(i,j)
576  enddo ; enddo
577  do j=jsq,jeq ; do i=is,ie
578  pfv(i,j,k) = -(m(i,j+1,k) - m(i,j,k)) * g%IdyCv(i,j)
579  enddo ; enddo
580  enddo
581  endif ! use_EOS
582 
583  if (present(eta)) then
584  if (cs%tides) then
585  ! eta is the sea surface height relative to a time-invariant geoid, for
586  ! comparison with what is used for eta in btstep. See how e was calculated
587  ! about 200 lines above.
588  !$OMP parallel do default(shared)
589  do j=jsq,jeq+1 ; do i=isq,ieq+1
590  eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
591  enddo ; enddo
592  else
593  !$OMP parallel do default(shared)
594  do j=jsq,jeq+1 ; do i=isq,ieq+1
595  eta(i,j) = e(i,j,1)*gv%Z_to_H
596  enddo ; enddo
597  endif
598  endif
599 
600  if (cs%id_PFu_bc>0) call post_data(cs%id_PFu_bc, cs%PFu_bc, cs%diag)
601  if (cs%id_PFv_bc>0) call post_data(cs%id_PFv_bc, cs%PFv_bc, cs%diag)
602  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
603 
604 end subroutine pressureforce_mont_bouss
605 
606 !> Determines the partial derivative of the acceleration due
607 !! to pressure forces with the free surface height.
608 subroutine set_pbce_bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
609  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
610  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
611  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface height [Z ~> m].
612  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
613  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
614  real, intent(in) :: rho0 !< The "Boussinesq" ocean density [kg m-3].
615  real, intent(in) :: gfs_scale !< Ratio between gravity applied to top
616  !! interface and the gravitational acceleration of
617  !! the planet [nondim]. Usually this ratio is 1.
618  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
619  intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due
620  !! to free surface height anomalies
621  !! [L2 T-2 H-1 ~> m s-2].
622  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
623  optional, intent(in) :: rho_star !< The layer densities (maybe compressibility
624  !! compensated), times g/rho_0 [L2 Z-1 T-2 ~> m s-2].
625 
626  ! Local variables
627  real :: ihtot(szi_(g)) ! The inverse of the sum of the layer thicknesses [H-1 ~> m-1 or m2 kg-1].
628  real :: press(szi_(g)) ! Interface pressure [Pa].
629  real :: t_int(szi_(g)) ! Interface temperature [degC].
630  real :: s_int(szi_(g)) ! Interface salinity [ppt].
631  real :: dr_dt(szi_(g)) ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
632  real :: dr_ds(szi_(g)) ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
633  real :: rho_in_situ(szi_(g)) ! In-situ density at the top of a layer [R ~> kg m-3].
634  real :: g_rho0 ! A scaled version of g_Earth / Rho0 [L2 m3 Z-1 T-2 kg-1 ~> m4 s-2 kg-1]
635  real :: rho0xg ! g_Earth * Rho0 [kg s-2 m-1 Z-1 ~> kg s-2 m-2]
636  logical :: use_eos ! If true, density is calculated from T & S using
637  ! an equation of state.
638  real :: z_neglect ! A thickness that is so small it is usually lost
639  ! in roundoff and can be neglected [Z ~> m].
640  integer :: isq, ieq, jsq, jeq, nz, i, j, k
641 
642  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
643 
644  rho0xg = rho0*us%L_T_to_m_s**2 * gv%g_Earth
645  g_rho0 = gv%g_Earth / gv%Rho0
646  use_eos = associated(tv%eqn_of_state)
647  z_neglect = gv%H_subroundoff*gv%H_to_Z
648 
649  if (use_eos) then
650  if (present(rho_star)) then
651  !$OMP parallel do default(shared) private(Ihtot)
652  do j=jsq,jeq+1
653  do i=isq,ieq+1
654  ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
655  pbce(i,j,1) = gfs_scale * rho_star(i,j,1) * gv%H_to_Z
656  enddo
657  do k=2,nz ; do i=isq,ieq+1
658  pbce(i,j,k) = pbce(i,j,k-1) + (rho_star(i,j,k)-rho_star(i,j,k-1)) * &
659  ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
660  enddo ; enddo
661  enddo ! end of j loop
662  else
663  !$OMP parallel do default(shared) private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS)
664  do j=jsq,jeq+1
665  do i=isq,ieq+1
666  ihtot(i) = gv%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
667  press(i) = -rho0xg*e(i,j,1)
668  enddo
669  call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
670  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
671  do i=isq,ieq+1
672  pbce(i,j,1) = g_rho0*(gfs_scale * rho_in_situ(i)) * gv%H_to_Z
673  enddo
674  do k=2,nz
675  do i=isq,ieq+1
676  press(i) = -rho0xg*e(i,j,k)
677  t_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
678  s_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
679  enddo
680  call calculate_density_derivs(t_int, s_int, press, dr_dt, dr_ds, &
681  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
682  do i=isq,ieq+1
683  pbce(i,j,k) = pbce(i,j,k-1) + g_rho0 * &
684  ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i)) * &
685  (dr_dt(i)*(tv%T(i,j,k)-tv%T(i,j,k-1)) + &
686  dr_ds(i)*(tv%S(i,j,k)-tv%S(i,j,k-1)))
687  enddo
688  enddo
689  enddo ! end of j loop
690  endif
691  else ! not use_EOS
692  !$OMP parallel do default(shared) private(Ihtot)
693  do j=jsq,jeq+1
694  do i=isq,ieq+1
695  ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
696  pbce(i,j,1) = gv%g_prime(1) * gv%H_to_Z
697  enddo
698  do k=2,nz ; do i=isq,ieq+1
699  pbce(i,j,k) = pbce(i,j,k-1) + &
700  (gv%g_prime(k)*gv%H_to_Z) * ((e(i,j,k) - e(i,j,nz+1)) * ihtot(i))
701  enddo ; enddo
702  enddo ! end of j loop
703  endif ! use_EOS
704 
705 end subroutine set_pbce_bouss
706 
707 !> Determines the partial derivative of the acceleration due
708 !! to pressure forces with the column mass.
709 subroutine set_pbce_nonbouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
710  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
711  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
712  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: p !< Interface pressures [Pa].
713  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
714  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
715  real, intent(in) :: gfs_scale !< Ratio between gravity applied to top
716  !! interface and the gravitational acceleration of
717  !! the planet [nondim]. Usually this ratio is 1.
718  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: pbce !< The baroclinic pressure anomaly in each layer due
719  !! to free surface height anomalies
720  !! [L2 H-1 T-2 ~> m4 kg-1 s-2].
721  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: alpha_star !< The layer specific volumes
722  !! (maybe compressibility compensated) [R-1 ~> m3 kg-1].
723  ! Local variables
724  real, dimension(SZI_(G),SZJ_(G)) :: &
725  dpbce, & ! A barotropic correction to the pbce to enable the use of
726  ! a reduced gravity form of the equations [L2 H-1 T-2 ~> m4 kg-1 s-2].
727  c_htot ! dP_dH divided by the total ocean pressure [R L2 T-2 H-1 Pa-1 ~> m2 kg-1].
728  real :: t_int(szi_(g)) ! Interface temperature [degC].
729  real :: s_int(szi_(g)) ! Interface salinity [ppt].
730  real :: dr_dt(szi_(g)) ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
731  real :: dr_ds(szi_(g)) ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
732  real :: rho_in_situ(szi_(g)) ! In-situ density at an interface [R ~> kg m-3].
733  real :: alpha_lay(szk_(g)) ! The specific volume of each layer [R-1 ~> m3 kg-1].
734  real :: dalpha_int(szk_(g)+1) ! The change in specific volume across each interface [R-1 ~> m3 kg-1].
735  real :: dp_dh ! A factor that converts from thickness to pressure times other dimensional
736  ! conversion factors [R L2 T-2 H-1 ~> Pa m2 kg-1].
737  real :: dp_neglect ! A thickness that is so small it is usually lost
738  ! in roundoff and can be neglected [Pa].
739  logical :: use_eos ! If true, density is calculated from T & S using
740  ! an equation of state.
741  integer :: isq, ieq, jsq, jeq, nz, i, j, k
742 
743  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
744 
745  use_eos = associated(tv%eqn_of_state)
746 
747  dp_dh = gv%g_Earth * gv%H_to_RZ
748  dp_neglect = gv%H_to_Pa * gv%H_subroundoff
749 
750  if (use_eos) then
751  if (present(alpha_star)) then
752  !$OMP parallel do default(shared)
753  do j=jsq,jeq+1
754  do i=isq,ieq+1
755  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
756  pbce(i,j,nz) = dp_dh * alpha_star(i,j,nz)
757  enddo
758  do k=nz-1,1,-1 ; do i=isq,ieq+1
759  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1)) * c_htot(i,j)) * &
760  (alpha_star(i,j,k) - alpha_star(i,j,k+1))
761  enddo ; enddo
762  enddo
763  else
764  !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
765  do j=jsq,jeq+1
766  call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), &
767  rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
768  do i=isq,ieq+1
769  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
770  pbce(i,j,nz) = dp_dh / (rho_in_situ(i))
771  enddo
772  do k=nz-1,1,-1
773  do i=isq,ieq+1
774  t_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
775  s_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
776  enddo
777  call calculate_density(t_int, s_int, p(:,j,k+1), rho_in_situ, &
778  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
779  call calculate_density_derivs(t_int, s_int, p(:,j,k+1), dr_dt, dr_ds, &
780  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
781  do i=isq,ieq+1
782  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
783  ((dr_dt(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &
784  dr_ds(i)*(tv%S(i,j,k+1)-tv%S(i,j,k))) / (rho_in_situ(i)**2))
785  enddo
786  enddo
787  enddo
788  endif
789  else ! not use_EOS
790 
791  do k=1,nz ; alpha_lay(k) = 1.0 / (gv%Rlay(k)) ; enddo
792  do k=2,nz ; dalpha_int(k) = alpha_lay(k-1) - alpha_lay(k) ; enddo
793 
794  !$OMP parallel do default(shared)
795  do j=jsq,jeq+1
796  do i=isq,ieq+1
797  c_htot(i,j) = dp_dh / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
798  pbce(i,j,nz) = dp_dh * alpha_lay(nz)
799  enddo
800  do k=nz-1,1,-1 ; do i=isq,ieq+1
801  pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,k+1)-p(i,j,1))*c_htot(i,j)) * &
802  dalpha_int(k+1)
803  enddo ; enddo
804  enddo
805  endif ! use_EOS
806 
807  if (gfs_scale < 1.0) then
808  ! Adjust the Montgomery potential to make this a reduced gravity model.
809  !$OMP parallel do default(shared)
810  do j=jsq,jeq+1 ; do i=isq,ieq+1
811  dpbce(i,j) = (gfs_scale - 1.0) * pbce(i,j,1)
812  enddo ; enddo
813  !$OMP parallel do default(shared)
814  do k=1,nz ; do j=jsq,jeq+1 ; do i=isq,ieq+1
815  pbce(i,j,k) = pbce(i,j,k) + dpbce(i,j)
816  enddo ; enddo ; enddo
817  endif
818 
819 end subroutine set_pbce_nonbouss
820 
821 !> Initialize the Montgomery-potential form of PGF control structure
822 subroutine pressureforce_mont_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
823  type(time_type), target, intent(in) :: time !< Current model time
824  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
825  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
826  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
827  type(param_file_type), intent(in) :: param_file !< Parameter file handles
828  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
829  type(pressureforce_mont_cs), pointer :: cs !< Montgomery PGF control structure
830  type(tidal_forcing_cs), optional, pointer :: tides_csp !< Tides control structure
831 
832  ! Local variables
833  logical :: use_temperature, use_eos
834  ! This include declares and sets the variable "version".
835 # include "version_variable.h"
836  character(len=40) :: mdl ! This module's name.
837 
838  if (associated(cs)) then
839  call mom_error(warning, "PressureForce_init called with an associated "// &
840  "control structure.")
841  return
842  else ; allocate(cs) ; endif
843 
844  cs%diag => diag ; cs%Time => time
845  if (present(tides_csp)) then
846  if (associated(tides_csp)) cs%tides_CSp => tides_csp
847  endif
848 
849  mdl = "MOM_PressureForce_Mont"
850  call log_version(param_file, mdl, version, "")
851  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
852  "The mean ocean density used with BOUSSINESQ true to "//&
853  "calculate accelerations and the mass for conservation "//&
854  "properties, or with BOUSSINSEQ false to convert some "//&
855  "parameters from vertical units of m to kg m-2.", &
856  units="kg m-3", default=1035.0)
857  call get_param(param_file, mdl, "TIDES", cs%tides, &
858  "If true, apply tidal momentum forcing.", default=.false.)
859  call get_param(param_file, mdl, "USE_EOS", use_eos, default=.true., &
860  do_not_log=.true.) ! Input for diagnostic use only.
861 
862  if (use_eos) then
863  cs%id_PFu_bc = register_diag_field('ocean_model', 'PFu_bc', diag%axesCuL, time, &
864  'Density Gradient Zonal Pressure Force Accel.', "meter second-2", conversion=us%L_T2_to_m_s2)
865  cs%id_PFv_bc = register_diag_field('ocean_model', 'PFv_bc', diag%axesCvL, time, &
866  'Density Gradient Meridional Pressure Force Accel.', "meter second-2", conversion=us%L_T2_to_m_s2)
867  if (cs%id_PFu_bc > 0) then
868  call safe_alloc_ptr(cs%PFu_bc,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
869  cs%PFu_bc(:,:,:) = 0.0
870  endif
871  if (cs%id_PFv_bc > 0) then
872  call safe_alloc_ptr(cs%PFv_bc,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
873  cs%PFv_bc(:,:,:) = 0.0
874  endif
875  endif
876 
877  if (cs%tides) then
878  cs%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
879  time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=us%Z_to_m)
880  endif
881 
882  cs%GFS_scale = 1.0
883  if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
884 
885  call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale)
886 
887 end subroutine pressureforce_mont_init
888 
889 !> Deallocates the Montgomery-potential form of PGF control structure
890 subroutine pressureforce_mont_end(CS)
891  type(pressureforce_mont_cs), pointer :: cs !< Control structure for Montgomery potential PGF
892  if (associated(cs)) deallocate(cs)
893 end subroutine pressureforce_mont_end
894 
895 !>\namespace mom_pressureforce_mont
896 !!
897 !! Provides the Boussunesq and non-Boussinesq forms of the horizontal
898 !! accelerations due to pressure gradients using the Montgomery potential. A
899 !! second-order accurate, centered scheme is used. If a split time stepping
900 !! scheme is used, the vertical decomposition into barotropic and baroclinic
901 !! contributions described by Hallberg (J Comp Phys 1997) is used. With a
902 !! nonlinear equation of state, compressibility is added along the lines proposed
903 !! by Sun et al. (JPO 1999), but with compressibility coefficients based on a fit
904 !! to a user-provided reference profile.
905 
906 end module mom_pressureforce_mont
mom_eos::int_specific_vol_dp
subroutine, public int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp)
Calls the appropriate subroutine to alculate analytical and nearly-analytical integrals in pressure a...
Definition: MOM_EOS.F90:665
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_eos::query_compressible
logical function, public query_compressible(EOS)
Returns true if the equation of state is compressible (i.e. has pressure dependence)
Definition: MOM_EOS.F90:796
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_pressureforce_mont::pressureforce_mont_nonbouss
subroutine, public pressureforce_mont_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
Non-Boussinesq Montgomery-potential form of pressure gradient.
Definition: MOM_PressureForce_Montgomery.F90:65
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_tidal_forcing
Tidal contributions to geopotential.
Definition: MOM_tidal_forcing.F90:2
mom_tidal_forcing::calc_tidal_forcing
subroutine, public calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to_Z)
This subroutine calculates the geopotential anomalies that drive the tides, including self-attraction...
Definition: MOM_tidal_forcing.F90:400
mom_tidal_forcing::tidal_forcing_cs
The control structure for the MOM_tidal_forcing module.
Definition: MOM_tidal_forcing.F90:26
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_pressureforce_mont::pressureforce_mont_bouss
subroutine, public pressureforce_mont_bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, eta)
Boussinesq Montgomery-potential form of pressure gradient.
Definition: MOM_PressureForce_Montgomery.F90:366
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_pressureforce_mont::pressureforce_mont_init
subroutine, public pressureforce_mont_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
Initialize the Montgomery-potential form of PGF control structure.
Definition: MOM_PressureForce_Montgomery.F90:823
mom_pressureforce_mont
Provides the Montgomery potential form of pressure gradient.
Definition: MOM_PressureForce_Montgomery.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_pressureforce_mont::set_pbce_nonbouss
subroutine, public set_pbce_nonbouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
Determines the partial derivative of the acceleration due to pressure forces with the column mass.
Definition: MOM_PressureForce_Montgomery.F90:710
mom_pressureforce_mont::pressureforce_mont_cs
Control structure for the Montgomery potential form of pressure gradient.
Definition: MOM_PressureForce_Montgomery.F90:31
mom_pressureforce_mont::pressureforce_mont_end
subroutine, public pressureforce_mont_end(CS)
Deallocates the Montgomery-potential form of PGF control structure.
Definition: MOM_PressureForce_Montgomery.F90:891
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_pressureforce_mont::set_pbce_bouss
subroutine, public set_pbce_bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
Determines the partial derivative of the acceleration due to pressure forces with the free surface he...
Definition: MOM_PressureForce_Montgomery.F90:609
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60