MOM_PressureForce_Montgomery.F90

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