MOM6
MOM_PressureForce_analytic_FV.F90
Go to the documentation of this file.
1 !> Analytically integrated finite volume 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, fatal, warning, is_root_pe
10 use mom_grid, only : ocean_grid_type
22 
23 implicit none ; private
24 
25 #include <MOM_memory.h>
26 
29 
30 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34 
35 !> Finite volume pressure gradient control structure
36 type, public :: pressureforce_afv_cs ; private
37  logical :: tides !< If true, apply tidal momentum forcing.
38  real :: rho0 !< The density used in the Boussinesq
39  !! approximation [kg m-3].
40  real :: gfs_scale !< A scaling of the surface pressure gradients to
41  !! allow the use of a reduced gravity model [nondim].
42  type(time_type), pointer :: time !< A pointer to the ocean model's clock.
43  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
44  !! timing of diagnostic output.
45  logical :: usemasswghtinterp !< Use mass weighting in T/S interpolation
46  logical :: boundary_extrap !< Indicate whether high-order boundary
47  !! extrapolation should be used within boundary cells
48 
49  logical :: reconstruct !< If true, polynomial profiles of T & S will be
50  !! reconstructed and used in the integrals for the
51  !! finite volume pressure gradient calculation.
52  !! The default depends on whether regridding is being used.
53 
54  integer :: recon_scheme !< Order of the polynomial of the reconstruction of T & S
55  !! for the finite volume pressure gradient calculation.
56  !! By the default (1) is for a piecewise linear method
57 
58  integer :: id_e_tidal = -1 !< Diagnostic identifier
59  type(tidal_forcing_cs), pointer :: tides_csp => null() !< Tides control structure
60 end type pressureforce_afv_cs
61 
62 contains
63 
64 !> Thin interface between the model and the Boussinesq and non-Boussinesq
65 !! pressure force routines.
66 subroutine pressureforce_afv(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
67  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
68  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
69  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
70  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
71  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variables
72  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: pfu !< Zonal acceleration [L T-2 ~> m s-2]
73  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration [L T-2 ~> m s-2]
74  type(pressureforce_afv_cs), pointer :: cs !< Finite volume PGF control structure
75  type(ale_cs), pointer :: ale_csp !< ALE control structure
76  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean
77  !! or atmosphere-ocean interface [Pa].
78  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
79  !! anomaly in each layer due to eta anomalies
80  !! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
81  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
82  !! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
83  !! contributions or compressibility compensation.
84 
85  if (gv%Boussinesq) then
86  call pressureforce_afv_bouss(h, tv, pfu, pfv, g, gv, us, cs, ale_csp, p_atm, pbce, eta)
87  else
88  call pressureforce_afv_nonbouss(h, tv, pfu, pfv, g, gv, us, cs, ale_csp, p_atm, pbce, eta)
89  endif
90 
91 end subroutine pressureforce_afv
92 
93 !> \brief Non-Boussinesq analytically-integrated finite volume form of pressure gradient
94 !!
95 !! Determines the acceleration due to hydrostatic pressure forces, using
96 !! the analytic finite volume form of the Pressure gradient, and does not
97 !! make the Boussinesq approximation.
98 !!
99 !! To work, the following fields must be set outside of the usual (is:ie,js:je)
100 !! range before this subroutine is called:
101 !! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
102 subroutine pressureforce_afv_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
103  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
104  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
105  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
106  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> kg/m2]
107  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
108  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: pfu !< Zonal acceleration [L T-2 ~> m s-2]
109  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration [L T-2 ~> m s-2]
110  type(pressureforce_afv_cs), pointer :: cs !< Finite volume PGF control structure
111  type(ale_cs), pointer :: ale_csp !< ALE control structure
112  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean
113  !! or atmosphere-ocean interface [Pa].
114  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
115  !! anomaly in each layer due to eta anomalies
116  !! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
117  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
118  !! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
119  !! contributions or compressibility compensation.
120  ! Local variables
121  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: p ! Interface pressure [Pa].
122  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
123  t_tmp, & ! Temporary array of temperatures where layers that are lighter
124  ! than the mixed layer have the mixed layer's properties [degC].
125  s_tmp ! Temporary array of salinities where layers that are lighter
126  ! than the mixed layer have the mixed layer's properties [ppt].
127  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
128  s_t, & ! Top and bottom edge values for linear reconstructions
129  s_b, & ! of salinity within each layer [ppt].
130  t_t, & ! Top and bottom edge values for linear reconstructions
131  t_b ! of temperature within each layer [degC].
132  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
133  dza, & ! The change in geopotential anomaly between the top and bottom
134  ! of a layer [m2 s-2].
135  intp_dza ! The vertical integral in depth of the pressure anomaly less
136  ! the pressure anomaly at the top of the layer [Pa m2 s-2].
137  real, dimension(SZI_(G),SZJ_(G)) :: &
138  dp, & ! The (positive) change in pressure across a layer [Pa].
139  ssh, & ! The sea surface height anomaly, in depth units [Z ~> m].
140  e_tidal, & ! The bottom geopotential anomaly due to tidal forces from
141  ! astronomical sources and self-attraction and loading [Z ~> m].
142  dm, & ! The barotropic adjustment to the Montgomery potential to
143  ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
144  za ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the
145  ! interface atop a layer [m2 s-2].
146 
147  real, dimension(SZI_(G)) :: rho_cv_bl ! The coordinate potential density in the deepest variable
148  ! density near-surface layer [R ~> kg m-3].
149  real, dimension(SZIB_(G),SZJ_(G)) :: &
150  intx_za ! The zonal integral of the geopotential anomaly along the
151  ! interface below a layer, divided by the grid spacing [m2 s-2].
152  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
153  intx_dza ! The change in intx_za through a layer [m2 s-2].
154  real, dimension(SZI_(G),SZJB_(G)) :: &
155  inty_za ! The meridional integral of the geopotential anomaly along the
156  ! interface below a layer, divided by the grid spacing [m2 s-2].
157  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
158  inty_dza ! The change in inty_za through a layer [m2 s-2].
159  real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate
160  ! density, [Pa] (usually 2e7 Pa = 2000 dbar).
161 
162  real :: dp_neglect ! A thickness that is so small it is usually lost
163  ! in roundoff and can be neglected [Pa].
164  real :: g_earth_z ! A scaled version of g_Earth [m2 Z-1 s-2 ~> m s-2].
165  real :: i_gearth ! The inverse of g_Earth_z [s2 Z m-2 ~> s2 m-1]
166  real :: alpha_anom ! The in-situ specific volume, averaged over a
167  ! layer, less alpha_ref [m3 kg-1].
168  logical :: use_p_atm ! If true, use the atmospheric pressure.
169  logical :: use_ale ! If true, use an ALE pressure reconstruction.
170  logical :: use_eos ! If true, density is calculated from T & S using an
171  ! equation of state.
172  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
173 
174  real :: alpha_ref ! A reference specific volume [m3 kg-1], that is used
175  ! to reduce the impact of truncation errors.
176  real :: rho_in_situ(szi_(g)) ! The in situ density [kg m-3].
177  real :: pa_to_h ! A factor to convert from Pa to the thicknesss units (H).
178 ! real :: oneatm = 101325.0 ! 1 atm in [Pa] = [kg m-1 s-2]
179  real, parameter :: c1_6 = 1.0/6.0
180  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
181  integer :: i, j, k
182 
183  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
184  nkmb=gv%nk_rho_varies
185  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
186 
187  if (.not.associated(cs)) call mom_error(fatal, &
188  "MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.")
189 
190  use_p_atm = .false.
191  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
192  use_eos = associated(tv%eqn_of_state)
193  use_ale = .false.
194  if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
195 
196  dp_neglect = gv%H_to_Pa * gv%H_subroundoff
197  alpha_ref = 1.0/cs%Rho0
198  g_earth_z = us%L_T_to_m_s**2 * gv%g_Earth
199  i_gearth = 1.0 / g_earth_z
200 
201  if (use_p_atm) then
202  !$OMP parallel do default(shared)
203  do j=jsq,jeq+1 ; do i=isq,ieq+1
204  p(i,j,1) = p_atm(i,j)
205  enddo ; enddo
206  else
207  !$OMP parallel do default(shared)
208  do j=jsq,jeq+1 ; do i=isq,ieq+1
209  p(i,j,1) = 0.0 ! or oneatm
210  enddo ; enddo
211  endif
212  !$OMP parallel do default(shared)
213  do j=jsq,jeq+1 ; do k=2,nz+1 ; do i=isq,ieq+1
214  p(i,j,k) = p(i,j,k-1) + gv%H_to_Pa * h(i,j,k-1)
215  enddo ; enddo ; enddo
216 
217  if (use_eos) then
218  ! With a bulk mixed layer, replace the T & S of any layers that are
219  ! lighter than the the buffer layer with the properties of the buffer
220  ! layer. These layers will be massless anyway, and it avoids any
221  ! formal calculations with hydrostatically unstable profiles.
222  if (nkmb>0) then
223  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
224  tv_tmp%eqn_of_state => tv%eqn_of_state
225  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
226  !$OMP parallel do default(shared) private(Rho_cv_BL)
227  do j=jsq,jeq+1
228  do k=1,nkmb ; do i=isq,ieq+1
229  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
230  enddo ; enddo
231  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
232  rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
233  do k=nkmb+1,nz ; do i=isq,ieq+1
234  if (gv%Rlay(k) < rho_cv_bl(i)) then
235  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
236  else
237  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
238  endif
239  enddo ; enddo
240  enddo
241  else
242  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
243  tv_tmp%eqn_of_state => tv%eqn_of_state
244  endif
245  endif
246 
247  ! If regridding is activated, do a linear reconstruction of salinity
248  ! and temperature across each layer. The subscripts 't' and 'b' refer
249  ! to top and bottom values within each layer (these are the only degrees
250  ! of freedeom needed to know the linear profile).
251  if ( use_ale ) then
252  if ( cs%Recon_Scheme == 1 ) then
253  call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
254  elseif ( cs%Recon_Scheme == 2) then
255  call pressure_gradient_ppm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
256  endif
257  endif
258 
259  !$OMP parallel do default(shared) private(alpha_anom,dp)
260  do k=1,nz
261  ! Calculate 4 integrals through the layer that are required in the
262  ! subsequent calculation.
263  if (use_eos) then
264  if ( use_ale ) then
265  if ( cs%Recon_Scheme == 1 ) then
266  call int_spec_vol_dp_generic_plm( t_t(:,:,k), t_b(:,:,k), &
267  s_t(:,:,k), s_b(:,:,k), p(:,:,k), p(:,:,k+1), &
268  alpha_ref, dp_neglect, p(:,:,nz+1), g%HI, &
269  tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
270  intx_dza(:,:,k), inty_dza(:,:,k), &
271  usemasswghtinterp = cs%useMassWghtInterp)
272  i=k
273  elseif ( cs%Recon_Scheme == 2 ) then
274  call mom_error(fatal, "PressureForce_AFV_nonBouss: "//&
275  "int_spec_vol_dp_generic_ppm does not exist yet.")
276  ! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
277  ! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), &
278  ! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
279  ! intx_dza(:,:,k), inty_dza(:,:,k))
280  endif
281  else
282  call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,k), &
283  p(:,:,k+1), alpha_ref, g%HI, tv%eqn_of_state, &
284  dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
285  inty_dza(:,:,k), bathyp=p(:,:,nz+1), dp_tiny=dp_neglect, &
286  usemasswghtinterp = cs%useMassWghtInterp)
287  endif
288  else
289  alpha_anom = 1.0/(us%R_to_kg_m3*gv%Rlay(k)) - alpha_ref
290  do j=jsq,jeq+1 ; do i=isq,ieq+1
291  dp(i,j) = gv%H_to_Pa * h(i,j,k)
292  dza(i,j,k) = alpha_anom * dp(i,j)
293  intp_dza(i,j,k) = 0.5 * alpha_anom * dp(i,j)**2
294  enddo ; enddo
295  do j=js,je ; do i=isq,ieq
296  intx_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i+1,j))
297  enddo ; enddo
298  do j=jsq,jeq ; do i=is,ie
299  inty_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i,j+1))
300  enddo ; enddo
301  endif
302  enddo
303 
304  ! The bottom geopotential anomaly is calculated first so that the increments
305  ! to the geopotential anomalies can be reused. Alternately, the surface
306  ! geopotential could be calculated directly with separate calls to
307  ! int_specific_vol_dp with alpha_ref=0, and the anomalies used going
308  ! downward, which would relieve the need for dza, intp_dza, intx_dza, and
309  ! inty_dza to be 3-D arrays.
310 
311  ! Sum vertically to determine the surface geopotential anomaly.
312  !$OMP parallel do default(shared)
313  do j=jsq,jeq+1
314  do i=isq,ieq+1
315  za(i,j) = alpha_ref*p(i,j,nz+1) - g_earth_z*g%bathyT(i,j)
316  enddo
317  do k=nz,1,-1 ; do i=isq,ieq+1
318  za(i,j) = za(i,j) + dza(i,j,k)
319  enddo ; enddo
320  enddo
321 
322  if (cs%tides) then
323  ! Find and add the tidal geopotential anomaly.
324  !$OMP parallel do default(shared)
325  do j=jsq,jeq+1 ; do i=isq,ieq+1
326  ssh(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * i_gearth
327  enddo ; enddo
328  call calc_tidal_forcing(cs%Time, ssh, e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
329  !$OMP parallel do default(shared)
330  do j=jsq,jeq+1 ; do i=isq,ieq+1
331  za(i,j) = za(i,j) - g_earth_z * e_tidal(i,j)
332  enddo ; enddo
333  endif
334 
335  if (cs%GFS_scale < 1.0) then
336  ! Adjust the Montgomery potential to make this a reduced gravity model.
337  if (use_eos) then
338  !$OMP parallel do default(shared) private(rho_in_situ)
339  do j=jsq,jeq+1
340  call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), &
341  rho_in_situ, isq, ieq-isq+2, tv%eqn_of_state)
342 
343  do i=isq,ieq+1
344  dm(i,j) = (cs%GFS_scale - 1.0) * us%m_s_to_L_T**2 * &
345  (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
346  enddo
347  enddo
348  else
349  !$OMP parallel do default(shared)
350  do j=jsq,jeq+1 ; do i=isq,ieq+1
351  dm(i,j) = (cs%GFS_scale - 1.0) * us%m_s_to_L_T**2 * &
352  (p(i,j,1)*(1.0/(us%R_to_kg_m3*gv%Rlay(1)) - alpha_ref) + za(i,j))
353  enddo ; enddo
354  endif
355 ! else
356 ! do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; dM(i,j) = 0.0 ; enddo ; enddo
357  endif
358 
359  ! This order of integrating upward and then downward again is necessary with
360  ! a nonlinear equation of state, so that the surface geopotentials will go
361  ! linearly between the values at thickness points, but the bottom
362  ! geopotentials will not now be linear at the sub-grid-scale. Doing this
363  ! ensures no motion with flat isopycnals, even with a nonlinear equation of state.
364  !$OMP parallel do default(shared)
365  do j=js,je ; do i=isq,ieq
366  intx_za(i,j) = 0.5*(za(i,j) + za(i+1,j))
367  enddo ; enddo
368  !$OMP parallel do default(shared)
369  do j=jsq,jeq ; do i=is,ie
370  inty_za(i,j) = 0.5*(za(i,j) + za(i,j+1))
371  enddo ; enddo
372  do k=1,nz
373  ! These expressions for the acceleration have been carefully checked in
374  ! a set of idealized cases, and should be bug-free.
375  !$OMP parallel do default(shared)
376  do j=jsq,jeq+1 ; do i=isq,ieq+1
377  dp(i,j) = gv%H_to_Pa*h(i,j,k)
378  za(i,j) = za(i,j) - dza(i,j,k)
379  enddo ; enddo
380  !$OMP parallel do default(shared)
381  do j=js,je ; do i=isq,ieq
382  intx_za(i,j) = intx_za(i,j) - intx_dza(i,j,k)
383  pfu(i,j,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
384  (za(i+1,j)*dp(i+1,j) + intp_dza(i+1,j,k))) + &
385  ((dp(i+1,j) - dp(i,j)) * intx_za(i,j) - &
386  (p(i+1,j,k) - p(i,j,k)) * intx_dza(i,j,k))) * &
387  (us%m_s_to_L_T**2 * 2.0*g%IdxCu(i,j) / &
388  ((dp(i,j) + dp(i+1,j)) + dp_neglect))
389  enddo ; enddo
390  !$OMP parallel do default(shared)
391  do j=jsq,jeq ; do i=is,ie
392  inty_za(i,j) = inty_za(i,j) - inty_dza(i,j,k)
393  pfv(i,j,k) = (((za(i,j)*dp(i,j) + intp_dza(i,j,k)) - &
394  (za(i,j+1)*dp(i,j+1) + intp_dza(i,j+1,k))) + &
395  ((dp(i,j+1) - dp(i,j)) * inty_za(i,j) - &
396  (p(i,j+1,k) - p(i,j,k)) * inty_dza(i,j,k))) * &
397  (us%m_s_to_L_T**2 * 2.0*g%IdyCv(i,j) / &
398  ((dp(i,j) + dp(i,j+1)) + dp_neglect))
399  enddo ; enddo
400 
401  if (cs%GFS_scale < 1.0) then
402  ! Adjust the Montgomery potential to make this a reduced gravity model.
403  !$OMP parallel do default(shared)
404  do j=js,je ; do i=isq,ieq
405  pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
406  enddo ; enddo
407  !$OMP parallel do default(shared)
408  do j=jsq,jeq ; do i=is,ie
409  pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
410  enddo ; enddo
411  endif
412  enddo
413 
414  if (present(pbce)) then
415  call set_pbce_nonbouss(p, tv_tmp, g, gv, us, cs%GFS_scale, pbce)
416  endif
417 
418  if (present(eta)) then
419  pa_to_h = 1.0 / gv%H_to_Pa
420  if (use_p_atm) then
421  !$OMP parallel do default(shared)
422  do j=jsq,jeq+1 ; do i=isq,ieq+1
423  eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*pa_to_h ! eta has the same units as h.
424  enddo ; enddo
425  else
426  !$OMP parallel do default(shared)
427  do j=jsq,jeq+1 ; do i=isq,ieq+1
428  eta(i,j) = p(i,j,nz+1)*pa_to_h ! eta has the same units as h.
429  enddo ; enddo
430  endif
431  endif
432 
433  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
434 
435 end subroutine pressureforce_afv_nonbouss
436 
437 !> \brief Boussinesq analytically-integrated finite volume form of pressure gradient
438 !!
439 !! Determines the acceleration due to hydrostatic pressure forces, using
440 !! the finite volume form of the terms and analytic integrals in depth.
441 !!
442 !! To work, the following fields must be set outside of the usual (is:ie,js:je)
443 !! range before this subroutine is called:
444 !! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
445 subroutine pressureforce_afv_bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
446  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
447  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
448  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
449  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m]
450  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
451  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: pfu !< Zonal acceleration [L T-2 ~> m s-2]
452  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: pfv !< Meridional acceleration [L T-2 ~> m s-2]
453  type(pressureforce_afv_cs), pointer :: cs !< Finite volume PGF control structure
454  type(ale_cs), pointer :: ale_csp !< ALE control structure
455  real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean
456  !! or atmosphere-ocean interface [Pa].
457  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
458  !! anomaly in each layer due to eta anomalies
459  !! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
460  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
461  !! calculate PFu and PFv [H ~> m or kg m-2], with any
462  !! tidal contributions or compressibility compensation.
463  ! Local variables
464  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: e ! Interface height in depth units [Z ~> m].
465  real, dimension(SZI_(G),SZJ_(G)) :: &
466  e_tidal, & ! The bottom geopotential anomaly due to tidal forces from
467  ! astronomical sources and self-attraction and loading [Z ~> m].
468  dm ! The barotropic adjustment to the Montgomery potential to
469  ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
470  real, dimension(SZI_(G)) :: &
471  rho_cv_bl ! The coordinate potential density in the deepest variable
472  ! density near-surface layer [R ~> kg m-3].
473  real, dimension(SZI_(G),SZJ_(G)) :: &
474  dz_geo, & ! The change in geopotential thickness through a layer times some dimensional
475  ! rescaling factors [kg m-1 R-1 s-2 ~> m2 s-2].
476  pa, & ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the
477  ! the interface atop a layer [Pa].
478  dpa, & ! The change in pressure anomaly between the top and bottom
479  ! of a layer [Pa].
480  intz_dpa ! The vertical integral in depth of the pressure anomaly less the
481  ! pressure anomaly at the top of the layer [H Pa ~> m Pa or kg m-2 Pa].
482  real, dimension(SZIB_(G),SZJ_(G)) :: &
483  intx_pa, & ! The zonal integral of the pressure anomaly along the interface
484  ! atop a layer, divided by the grid spacing [Pa].
485  intx_dpa ! The change in intx_pa through a layer [Pa].
486  real, dimension(SZI_(G),SZJB_(G)) :: &
487  inty_pa, & ! The meridional integral of the pressure anomaly along the
488  ! interface atop a layer, divided by the grid spacing [Pa].
489  inty_dpa ! The change in inty_pa through a layer [Pa].
490 
491  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
492  t_tmp, & ! Temporary array of temperatures where layers that are lighter
493  ! than the mixed layer have the mixed layer's properties [degC].
494  s_tmp ! Temporary array of salinities where layers that are lighter
495  ! than the mixed layer have the mixed layer's properties [ppt].
496  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
497  s_t, s_b, t_t, t_b ! Top and bottom edge values for linear reconstructions
498  ! of salinity and temperature within each layer.
499  real :: rho_in_situ(szi_(g)) ! The in situ density [R ~> kg m-3].
500  real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate
501  ! density, [Pa] (usually 2e7 Pa = 2000 dbar).
502  real :: p0(szi_(g)) ! An array of zeros to use for pressure [Pa].
503  real :: h_neglect ! A thickness that is so small it is usually lost
504  ! in roundoff and can be neglected [H ~> m].
505  real :: g_earth_mks_z ! A scaled version of g_Earth [m2 Z-1 s-2 ~> m s-2].
506  real :: g_earth_z_geo ! Another scaled version of g_Earth [R m5 kg-1 Z-1 s-2 ~> m s-2].
507  real :: i_rho0 ! 1/Rho0 times unit scaling factors [L2 m kg-1 s2 T-2 ~> m3 kg-1].
508  real :: g_rho0 ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
509  real :: rho_ref ! The reference density [R ~> kg m-3].
510  real :: rho_ref_mks ! The reference density in mks units [kg m-3].
511  real :: dz_neglect ! A minimal thickness [Z ~> m], like e.
512  logical :: use_p_atm ! If true, use the atmospheric pressure.
513  logical :: use_ale ! If true, use an ALE pressure reconstruction.
514  logical :: use_eos ! If true, density is calculated from T & S using an equation of state.
515  type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
516 
517  real, parameter :: c1_6 = 1.0/6.0
518  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
519  integer :: i, j, k
520 
521  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
522  nkmb=gv%nk_rho_varies
523  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
524 
525  if (.not.associated(cs)) call mom_error(fatal, &
526  "MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.")
527 
528  use_p_atm = .false.
529  if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
530  use_eos = associated(tv%eqn_of_state)
531  do i=isq,ieq+1 ; p0(i) = 0.0 ; enddo
532  use_ale = .false.
533  if (associated(ale_csp)) use_ale = cs%reconstruct .and. use_eos
534 
535  h_neglect = gv%H_subroundoff
536  dz_neglect = gv%H_subroundoff * gv%H_to_Z
537  i_rho0 = us%m_s_to_L_T**2 / (us%R_to_kg_m3*gv%Rho0)
538  g_earth_mks_z = us%L_T_to_m_s**2 * gv%g_Earth
539  g_earth_z_geo = us%R_to_kg_m3*us%L_T_to_m_s**2 * gv%g_Earth
540  g_rho0 = gv%g_Earth / gv%Rho0
541  rho_ref_mks = cs%Rho0
542  rho_ref = rho_ref_mks*us%kg_m3_to_R
543 
544  if (cs%tides) then
545  ! Determine the surface height anomaly for calculating self attraction
546  ! and loading. This should really be based on bottom pressure anomalies,
547  ! but that is not yet implemented, and the current form is correct for
548  ! barotropic tides.
549  !$OMP parallel do default(shared)
550  do j=jsq,jeq+1
551  do i=isq,ieq+1
552  e(i,j,1) = -g%bathyT(i,j)
553  enddo
554  do k=1,nz ; do i=isq,ieq+1
555  e(i,j,1) = e(i,j,1) + h(i,j,k)*gv%H_to_Z
556  enddo ; enddo
557  enddo
558  call calc_tidal_forcing(cs%Time, e(:,:,1), e_tidal, g, cs%tides_CSp, m_to_z=us%m_to_Z)
559  endif
560 
561 ! Here layer interface heights, e, are calculated.
562  if (cs%tides) then
563  !$OMP parallel do default(shared)
564  do j=jsq,jeq+1 ; do i=isq,ieq+1
565  e(i,j,nz+1) = -(g%bathyT(i,j) + e_tidal(i,j))
566  enddo ; enddo
567  else
568  !$OMP parallel do default(shared)
569  do j=jsq,jeq+1 ; do i=isq,ieq+1
570  e(i,j,nz+1) = -g%bathyT(i,j)
571  enddo ; enddo
572  endif
573  !$OMP parallel do default(shared)
574  do j=jsq,jeq+1; do k=nz,1,-1 ; do i=isq,ieq+1
575  e(i,j,k) = e(i,j,k+1) + h(i,j,k)*gv%H_to_Z
576  enddo ; enddo ; enddo
577 
578  if (use_eos) then
579 ! With a bulk mixed layer, replace the T & S of any layers that are
580 ! lighter than the the buffer layer with the properties of the buffer
581 ! layer. These layers will be massless anyway, and it avoids any
582 ! formal calculations with hydrostatically unstable profiles.
583 
584  if (nkmb>0) then
585  tv_tmp%T => t_tmp ; tv_tmp%S => s_tmp
586  tv_tmp%eqn_of_state => tv%eqn_of_state
587 
588  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
589  !$OMP parallel do default(shared) private(Rho_cv_BL)
590  do j=jsq,jeq+1
591  do k=1,nkmb ; do i=isq,ieq+1
592  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
593  enddo ; enddo
594  call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
595  rho_cv_bl(:), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
596 
597  do k=nkmb+1,nz ; do i=isq,ieq+1
598  if (gv%Rlay(k) < rho_cv_bl(i)) then
599  tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
600  else
601  tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
602  endif
603  enddo ; enddo
604  enddo
605  else
606  tv_tmp%T => tv%T ; tv_tmp%S => tv%S
607  tv_tmp%eqn_of_state => tv%eqn_of_state
608  endif
609  endif
610 
611  if (cs%GFS_scale < 1.0) then
612  ! Adjust the Montgomery potential to make this a reduced gravity model.
613  if (use_eos) then
614  !$OMP parallel do default(shared)
615  do j=jsq,jeq+1
616  if (use_p_atm) then
617  call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, &
618  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
619  else
620  call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, &
621  isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
622  endif
623  do i=isq,ieq+1
624  dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * rho_in_situ(i)) * e(i,j,1)
625  enddo
626  enddo
627  else
628  !$OMP parallel do default(shared)
629  do j=jsq,jeq+1 ; do i=isq,ieq+1
630  dm(i,j) = (cs%GFS_scale - 1.0) * (g_rho0 * gv%Rlay(1)) * e(i,j,1)
631  enddo ; enddo
632  endif
633  endif
634  ! I have checked that rho_0 drops out and that the 1-layer case is right. RWH.
635 
636  ! If regridding is activated, do a linear reconstruction of salinity
637  ! and temperature across each layer. The subscripts 't' and 'b' refer
638  ! to top and bottom values within each layer (these are the only degrees
639  ! of freedeom needed to know the linear profile).
640  if ( use_ale ) then
641  if ( cs%Recon_Scheme == 1 ) then
642  call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
643  elseif ( cs%Recon_Scheme == 2 ) then
644  call pressure_gradient_ppm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, cs%boundary_extrap)
645  endif
646  endif
647 
648  ! Set the surface boundary conditions on pressure anomaly and its horizontal
649  ! integrals, assuming that the surface pressure anomaly varies linearly
650  ! in x and y.
651  if (use_p_atm) then
652  !$OMP parallel do default(shared)
653  do j=jsq,jeq+1 ; do i=isq,ieq+1
654  pa(i,j) = (rho_ref*g_earth_z_geo)*e(i,j,1) + p_atm(i,j)
655  enddo ; enddo
656  else
657  !$OMP parallel do default(shared)
658  do j=jsq,jeq+1 ; do i=isq,ieq+1
659  pa(i,j) = (rho_ref*g_earth_z_geo)*e(i,j,1)
660  enddo ; enddo
661  endif
662  !$OMP parallel do default(shared)
663  do j=js,je ; do i=isq,ieq
664  intx_pa(i,j) = 0.5*(pa(i,j) + pa(i+1,j))
665  enddo ; enddo
666  !$OMP parallel do default(shared)
667  do j=jsq,jeq ; do i=is,ie
668  inty_pa(i,j) = 0.5*(pa(i,j) + pa(i,j+1))
669  enddo ; enddo
670 
671  do k=1,nz
672  ! Calculate 4 integrals through the layer that are required in the
673  ! subsequent calculation.
674 
675  if (use_eos) then
676  ! The following routine computes the integrals that are needed to
677  ! calculate the pressure gradient force. Linear profiles for T and S are
678  ! assumed when regridding is activated. Otherwise, the previous version
679  ! is used, whereby densities within each layer are constant no matter
680  ! where the layers are located.
681  if ( use_ale ) then
682  if ( cs%Recon_Scheme == 1 ) then
683  call int_density_dz_generic_plm( t_t(:,:,k), t_b(:,:,k), &
684  s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
685  rho_ref_mks, cs%Rho0, g_earth_mks_z, &
686  dz_neglect, g%bathyT, g%HI, g%HI, &
687  tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, &
688  usemasswghtinterp = cs%useMassWghtInterp)
689  elseif ( cs%Recon_Scheme == 2 ) then
690  call int_density_dz_generic_ppm( tv%T(:,:,k), t_t(:,:,k), t_b(:,:,k), &
691  tv%S(:,:,k), s_t(:,:,k), s_b(:,:,k), e(:,:,k), e(:,:,k+1), &
692  rho_ref_mks, cs%Rho0, g_earth_mks_z, &
693  g%HI, g%HI, tv%eqn_of_state, dpa, intz_dpa, &
694  intx_dpa, inty_dpa)
695  endif
696  else
697  call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,k), e(:,:,k+1), &
698  rho_ref_mks, cs%Rho0, g_earth_mks_z, g%HI, g%HI, tv%eqn_of_state, &
699  dpa, intz_dpa, intx_dpa, inty_dpa, &
700  g%bathyT, dz_neglect, cs%useMassWghtInterp)
701  endif
702  !$OMP parallel do default(shared)
703  do j=jsq,jeq+1 ; do i=isq,ieq+1
704  intz_dpa(i,j) = intz_dpa(i,j)*gv%Z_to_H
705  enddo ; enddo
706  else
707  !$OMP parallel do default(shared)
708  do j=jsq,jeq+1 ; do i=isq,ieq+1
709  dz_geo(i,j) = g_earth_z_geo * gv%H_to_Z*h(i,j,k)
710  dpa(i,j) = (gv%Rlay(k) - rho_ref) * dz_geo(i,j)
711  intz_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * dz_geo(i,j)*h(i,j,k)
712  enddo ; enddo
713  !$OMP parallel do default(shared)
714  do j=js,je ; do i=isq,ieq
715  intx_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i+1,j))
716  enddo ; enddo
717  !$OMP parallel do default(shared)
718  do j=jsq,jeq ; do i=is,ie
719  inty_dpa(i,j) = 0.5*(gv%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i,j+1))
720  enddo ; enddo
721  endif
722 
723  ! Compute pressure gradient in x direction
724  !$OMP parallel do default(shared)
725  do j=js,je ; do i=isq,ieq
726  pfu(i,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
727  (pa(i+1,j)*h(i+1,j,k) + intz_dpa(i+1,j))) + &
728  ((h(i+1,j,k) - h(i,j,k)) * intx_pa(i,j) - &
729  (e(i+1,j,k+1) - e(i,j,k+1)) * intx_dpa(i,j) * gv%Z_to_H)) * &
730  ((2.0*i_rho0*g%IdxCu(i,j)) / &
731  ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
732  intx_pa(i,j) = intx_pa(i,j) + intx_dpa(i,j)
733  enddo ; enddo
734  ! Compute pressure gradient in y direction
735  !$OMP parallel do default(shared)
736  do j=jsq,jeq ; do i=is,ie
737  pfv(i,j,k) = (((pa(i,j)*h(i,j,k) + intz_dpa(i,j)) - &
738  (pa(i,j+1)*h(i,j+1,k) + intz_dpa(i,j+1))) + &
739  ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,j) - &
740  (e(i,j+1,k+1) - e(i,j,k+1)) * inty_dpa(i,j) * gv%Z_to_H)) * &
741  ((2.0*i_rho0*g%IdyCv(i,j)) / &
742  ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
743  inty_pa(i,j) = inty_pa(i,j) + inty_dpa(i,j)
744  enddo ; enddo
745  !$OMP parallel do default(shared)
746  do j=jsq,jeq+1 ; do i=isq,ieq+1
747  pa(i,j) = pa(i,j) + dpa(i,j)
748  enddo ; enddo
749  enddo
750 
751  if (cs%GFS_scale < 1.0) then
752  do k=1,nz
753  !$OMP parallel do default(shared)
754  do j=js,je ; do i=isq,ieq
755  pfu(i,j,k) = pfu(i,j,k) - (dm(i+1,j) - dm(i,j)) * g%IdxCu(i,j)
756  enddo ; enddo
757  !$OMP parallel do default(shared)
758  do j=jsq,jeq ; do i=is,ie
759  pfv(i,j,k) = pfv(i,j,k) - (dm(i,j+1) - dm(i,j)) * g%IdyCv(i,j)
760  enddo ; enddo
761  enddo
762  endif
763 
764  if (present(pbce)) then
765  call set_pbce_bouss(e, tv_tmp, g, gv, us, cs%Rho0, cs%GFS_scale, pbce)
766  endif
767 
768  if (present(eta)) then
769  if (cs%tides) then
770  ! eta is the sea surface height relative to a time-invariant geoid, for
771  ! comparison with what is used for eta in btstep. See how e was calculated
772  ! about 200 lines above.
773  !$OMP parallel do default(shared)
774  do j=jsq,jeq+1 ; do i=isq,ieq+1
775  eta(i,j) = e(i,j,1)*gv%Z_to_H + e_tidal(i,j)*gv%Z_to_H
776  enddo ; enddo
777  else
778  !$OMP parallel do default(shared)
779  do j=jsq,jeq+1 ; do i=isq,ieq+1
780  eta(i,j) = e(i,j,1)*gv%Z_to_H
781  enddo ; enddo
782  endif
783  endif
784 
785  if (cs%id_e_tidal>0) call post_data(cs%id_e_tidal, e_tidal, cs%diag)
786 
787 end subroutine pressureforce_afv_bouss
788 
789 !> Initializes the finite volume pressure gradient control structure
790 subroutine pressureforce_afv_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
791  type(time_type), target, intent(in) :: time !< Current model time
792  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
793  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
794  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
795  type(param_file_type), intent(in) :: param_file !< Parameter file handles
796  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
797  type(pressureforce_afv_cs), pointer :: cs !< Finite volume PGF control structure
798  type(tidal_forcing_cs), optional, pointer :: tides_csp !< Tides control structure
799  ! This include declares and sets the variable "version".
800 # include "version_variable.h"
801  character(len=40) :: mdl ! This module's name.
802  logical :: use_ale
803 
804  if (associated(cs)) then
805  call mom_error(warning, "PressureForce_init called with an associated "// &
806  "control structure.")
807  return
808  else ; allocate(cs) ; endif
809 
810  cs%diag => diag ; cs%Time => time
811  if (present(tides_csp)) then
812  if (associated(tides_csp)) cs%tides_CSp => tides_csp
813  endif
814 
815  mdl = "MOM_PressureForce_AFV"
816  call log_version(param_file, mdl, version, "")
817  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
818  "The mean ocean density used with BOUSSINESQ true to "//&
819  "calculate accelerations and the mass for conservation "//&
820  "properties, or with BOUSSINSEQ false to convert some "//&
821  "parameters from vertical units of m to kg m-2.", &
822  units="kg m-3", default=1035.0)
823  call get_param(param_file, mdl, "TIDES", cs%tides, &
824  "If true, apply tidal momentum forcing.", default=.false.)
825  call get_param(param_file, "MOM", "USE_REGRIDDING", use_ale, &
826  "If True, use the ALE algorithm (regridding/remapping). "//&
827  "If False, use the layered isopycnal algorithm.", default=.false. )
828  call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", cs%useMassWghtInterp, &
829  "If true, use mass weighting when interpolating T/S for "//&
830  "integrals near the bathymetry in AFV pressure gradient "//&
831  "calculations.", default=.false.)
832  call get_param(param_file, mdl, "RECONSTRUCT_FOR_PRESSURE", cs%reconstruct, &
833  "If True, use vertical reconstruction of T & S within "//&
834  "the integrals of the FV pressure gradient calculation. "//&
835  "If False, use the constant-by-layer algorithm. "//&
836  "The default is set by USE_REGRIDDING.", &
837  default=use_ale )
838  call get_param(param_file, mdl, "PRESSURE_RECONSTRUCTION_SCHEME", cs%Recon_Scheme, &
839  "Order of vertical reconstruction of T/S to use in the "//&
840  "integrals within the FV pressure gradient calculation.\n"//&
841  " 0: PCM or no reconstruction.\n"//&
842  " 1: PLM reconstruction.\n"//&
843  " 2: PPM reconstruction.", default=1)
844  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION_PRESSURE", cs%boundary_extrap, &
845  "If true, the reconstruction of T & S for pressure in "//&
846  "boundary cells is extrapolated, rather than using PCM "//&
847  "in these cells. If true, the same order polynomial is "//&
848  "used as is used for the interior cells.", default=.true.)
849 
850  if (cs%tides) then
851  cs%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
852  time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=us%Z_to_m)
853  endif
854 
855  cs%GFS_scale = 1.0
856  if (gv%g_prime(1) /= gv%g_Earth) cs%GFS_scale = gv%g_prime(1) / gv%g_Earth
857 
858  call log_param(param_file, mdl, "GFS / G_EARTH", cs%GFS_scale)
859 
860 end subroutine pressureforce_afv_init
861 
862 !> Deallocates the finite volume pressure gradient control structure
863 subroutine pressureforce_afv_end(CS)
864  type(pressureforce_afv_cs), pointer :: cs !< Finite volume pressure control structure that
865  !! will be deallocated in this subroutine.
866  if (associated(cs)) deallocate(cs)
867 end subroutine pressureforce_afv_end
868 
869 !> \namespace mom_pressureforce_afv
870 !!
871 !! Provides the Boussinesq and non-Boussinesq forms of horizontal accelerations
872 !! due to pressure gradients using a 2nd-order analytically vertically integrated
873 !! finite volume form, as described by Adcroft et al., 2008.
874 !!
875 !! This form eliminates the thermobaric instabilities that had been a problem with
876 !! previous forms of the pressure gradient force calculation, as described by
877 !! Hallberg, 2005.
878 !!
879 !! Adcroft, A., R. Hallberg, and M. Harrison, 2008: A finite volume discretization
880 !! of the pressure gradient force using analytic integration. Ocean Modelling, 22,
881 !! 106-113. http://doi.org/10.1016/j.ocemod.2008.02.001
882 !!
883 !! Hallberg, 2005: A thermobaric instability of Lagrangian vertical coordinate
884 !! ocean models. Ocean Modelling, 8, 279-300.
885 !! http://dx.doi.org/10.1016/j.ocemod.2004.01.001
886 
887 end module mom_pressureforce_afv
mom_eos::int_density_dz_generic_plm
subroutine, public int_density_dz_generic_plm(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp)
Compute pressure gradient force integrals by quadrature for the case where T and S are linear profile...
Definition: MOM_EOS.F90:1171
mom_pressureforce_afv::pressureforce_afv
subroutine, public pressureforce_afv(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
Thin interface between the model and the Boussinesq and non-Boussinesq pressure force routines.
Definition: MOM_PressureForce_analytic_FV.F90:67
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_pressureforce_afv::pressureforce_afv_nonbouss
subroutine, public pressureforce_afv_nonbouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
Non-Boussinesq analytically-integrated finite volume form of pressure gradient.
Definition: MOM_PressureForce_analytic_FV.F90:103
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_eos::int_spec_vol_dp_generic
subroutine, public int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, useMassWghtInterp)
This subroutine calculates integrals of specific volume anomalies in pressure across layers,...
Definition: MOM_EOS.F90:2024
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_eos::int_density_dz_generic
subroutine, public int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
This subroutine calculates (by numerical quadrature) integrals of pressure anomalies across layers,...
Definition: MOM_EOS.F90:984
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_eos::int_density_dz
subroutine, public int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
This subroutine calculates analytical and nearly-analytical integrals of pressure anomalies across la...
Definition: MOM_EOS.F90:733
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_eos::int_density_dz_generic_ppm
subroutine, public int_density_dz_generic_ppm(T, T_t, T_b, S, S_t, S_b, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
Compute pressure gradient force integrals for the case where T and S are parabolic profiles.
Definition: MOM_EOS.F90:1585
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_eos::int_spec_vol_dp_generic_plm
subroutine, public int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, dP_neglect, bathyP, HI, EOS, dza, intp_dza, intx_dza, inty_dza, useMassWghtInterp)
This subroutine calculates integrals of specific volume anomalies in pressure across layers,...
Definition: MOM_EOS.F90:2214
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_pressureforce_afv::pressureforce_afv_cs
Finite volume pressure gradient control structure.
Definition: MOM_PressureForce_analytic_FV.F90:36
mom_ale
This module contains the main regridding routines.
Definition: MOM_ALE.F90:9
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_pressureforce_afv::pressureforce_afv_end
subroutine, public pressureforce_afv_end(CS)
Deallocates the finite volume pressure gradient control structure.
Definition: MOM_PressureForce_analytic_FV.F90:864
mom_ale::pressure_gradient_plm
subroutine, public pressure_gradient_plm(CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap)
Use plm reconstruction for pressure gradient (determine edge values) By using a PLM (limited piecewis...
Definition: MOM_ALE.F90:1012
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_ale::ale_cs
ALE control structure.
Definition: MOM_ALE.F90:65
mom_pressureforce_afv
Analytically integrated finite volume pressure gradient.
Definition: MOM_PressureForce_analytic_FV.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: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_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_pressureforce_afv::pressureforce_afv_init
subroutine, public pressureforce_afv_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
Initializes the finite volume pressure gradient control structure.
Definition: MOM_PressureForce_analytic_FV.F90:791
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_ale::pressure_gradient_ppm
subroutine, public pressure_gradient_ppm(CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap)
Use ppm reconstruction for pressure gradient (determine edge values) By using a PPM (limited piecewis...
Definition: MOM_ALE.F90:1086
mom_pressureforce_afv::pressureforce_afv_bouss
subroutine, public pressureforce_afv_bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
Boussinesq analytically-integrated finite volume form of pressure gradient.
Definition: MOM_PressureForce_analytic_FV.F90:446
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