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