MOM6
MOM_density_integrals.F90
1 !> Provides integrals of density
2 module mom_density_integrals
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_eos, only : eos_type
7 use mom_eos, only : eos_quadrature
8 use mom_eos, only : analytic_int_density_dz
9 use mom_eos, only : analytic_int_specific_vol_dp
10 use mom_eos, only : calculate_density
11 use mom_eos, only : calculate_spec_vol
12 use mom_eos, only : calculate_specific_vol_derivs
13 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg
14 use mom_file_parser, only : get_param, log_version, param_file_type
15 use mom_hor_index, only : hor_index_type
16 use mom_string_functions, only : uppercase
17 use mom_variables, only : thermo_var_ptrs
18 use mom_unit_scaling, only : unit_scale_type
19 use mom_verticalgrid, only : verticalgrid_type
20 
21 implicit none ; private
22 
23 #include <MOM_memory.h>
24 
25 public int_density_dz
26 public int_density_dz_generic_pcm
27 public int_density_dz_generic_plm
28 public int_density_dz_generic_ppm
29 public int_specific_vol_dp
30 public int_spec_vol_dp_generic_pcm
31 public int_spec_vol_dp_generic_plm
32 public find_depth_of_pressure_in_cell
33 
34 contains
35 
36 !> Calls the appropriate subroutine to calculate analytical and nearly-analytical
37 !! integrals in z across layers of pressure anomalies, which are
38 !! required for calculating the finite-volume form pressure accelerations in a
39 !! Boussinesq model.
40 subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, &
41  intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
42  type(hor_index_type), intent(in) :: hi !< Ocean horizontal index structures for the arrays
43  real, dimension(SZI_(HI),SZJ_(HI)), &
44  intent(in) :: t !< Potential temperature referenced to the surface [degC]
45  real, dimension(SZI_(HI),SZJ_(HI)), &
46  intent(in) :: s !< Salinity [ppt]
47  real, dimension(SZI_(HI),SZJ_(HI)), &
48  intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]
49  real, dimension(SZI_(HI),SZJ_(HI)), &
50  intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]
51  real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is
52  !! subtracted out to reduce the magnitude of each of the
53  !! integrals.
54  real, intent(in) :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used
55  !! to calculate the pressure (as p~=-z*rho_0*G_e)
56  !! used in the equation of state.
57  real, intent(in) :: g_e !< The Earth's gravitational acceleration
58  !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2]
59  type(eos_type), pointer :: eos !< Equation of state structure
60  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
61  real, dimension(SZI_(HI),SZJ_(HI)), &
62  intent(inout) :: dpa !< The change in the pressure anomaly
63  !! across the layer [R L2 T-2 ~> Pa] or [Pa]
64  real, dimension(SZI_(HI),SZJ_(HI)), &
65  optional, intent(inout) :: intz_dpa !< The integral through the thickness of the
66  !! layer of the pressure anomaly relative to the
67  !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]
68  real, dimension(SZIB_(HI),SZJ_(HI)), &
69  optional, intent(inout) :: intx_dpa !< The integral in x of the difference between
70  !! the pressure anomaly at the top and bottom of the
71  !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
72  real, dimension(SZI_(HI),SZJB_(HI)), &
73  optional, intent(inout) :: inty_dpa !< The integral in y of the difference between
74  !! the pressure anomaly at the top and bottom of the
75  !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
76  real, dimension(SZI_(HI),SZJ_(HI)), &
77  optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]
78  real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]
79  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
80  !! interpolate T/S for top and bottom integrals.
81 
82  if (eos_quadrature(eos)) then
83  call int_density_dz_generic_pcm(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, eos, us, dpa, &
84  intz_dpa, intx_dpa, inty_dpa, bathyt, dz_neglect, usemasswghtinterp)
85  else
86  call analytic_int_density_dz(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, eos, dpa, &
87  intz_dpa, intx_dpa, inty_dpa, bathyt, dz_neglect, usemasswghtinterp)
88  endif
89 
90 end subroutine int_density_dz
91 
92 
93 !> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which
94 !! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model.
95 subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
96  EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
97  bathyT, dz_neglect, useMassWghtInterp)
98  type(hor_index_type), intent(in) :: hi !< Horizontal index type for input variables.
99  real, dimension(SZI_(HI),SZJ_(HI)), &
100  intent(in) :: t !< Potential temperature of the layer [degC]
101  real, dimension(SZI_(HI),SZJ_(HI)), &
102  intent(in) :: s !< Salinity of the layer [ppt]
103  real, dimension(SZI_(HI),SZJ_(HI)), &
104  intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]
105  real, dimension(SZI_(HI),SZJ_(HI)), &
106  intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]
107  real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is
108  !! subtracted out to reduce the magnitude
109  !! of each of the integrals.
110  real, intent(in) :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used
111  !! to calculate the pressure (as p~=-z*rho_0*G_e)
112  !! used in the equation of state.
113  real, intent(in) :: g_e !< The Earth's gravitational acceleration
114  !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2]
115  type(eos_type), pointer :: eos !< Equation of state structure
116  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
117  real, dimension(SZI_(HI),SZJ_(HI)), &
118  intent(inout) :: dpa !< The change in the pressure anomaly
119  !! across the layer [R L2 T-2 ~> Pa] or [Pa]
120  real, dimension(SZI_(HI),SZJ_(HI)), &
121  optional, intent(inout) :: intz_dpa !< The integral through the thickness of the
122  !! layer of the pressure anomaly relative to the
123  !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]
124  real, dimension(SZIB_(HI),SZJ_(HI)), &
125  optional, intent(inout) :: intx_dpa !< The integral in x of the difference between
126  !! the pressure anomaly at the top and bottom of the
127  !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
128  real, dimension(SZI_(HI),SZJB_(HI)), &
129  optional, intent(inout) :: inty_dpa !< The integral in y of the difference between
130  !! the pressure anomaly at the top and bottom of the
131  !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
132  real, dimension(SZI_(HI),SZJ_(HI)), &
133  optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]
134  real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]
135  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
136  !! interpolate T/S for top and bottom integrals.
137  ! Local variables
138  real :: t5(5), s5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt]
139  real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa]
140  real :: r5(5) ! Densities at five quadrature points [R ~> kg m-3] or [kg m-3]
141  real :: rho_anom ! The depth averaged density anomaly [R ~> kg m-3] or [kg m-3]
142  real :: w_left, w_right ! Left and right weights [nondim]
143  real, parameter :: c1_90 = 1.0/90.0 ! Rational constants.
144  real :: gxrho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
145  real :: i_rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]
146  real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]
147  real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]
148  real :: dz ! The layer thickness [Z ~> m]
149  real :: hwght ! A pressure-thickness below topography [Z ~> m]
150  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m]
151  real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2]
152  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]
153  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]
154  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]
155  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]
156  real :: intz(5) ! The gravitational acceleration times the integrals of density
157  ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]
158  logical :: do_massweight ! Indicates whether to do mass weighting.
159  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n
160 
161  ! These array bounds work for the indexing convention of the input arrays, but
162  ! on the computational domain defined for the output arrays.
163  isq = hi%IscB ; ieq = hi%IecB
164  jsq = hi%JscB ; jeq = hi%JecB
165  is = hi%isc ; ie = hi%iec
166  js = hi%jsc ; je = hi%jec
167 
168  rho_scale = us%kg_m3_to_R
169  gxrho = us%RL2_T2_to_Pa * g_e * rho_0
170  rho_ref_mks = rho_ref * us%R_to_kg_m3
171  i_rho = 1.0 / rho_0
172 
173  do_massweight = .false.
174  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
175  do_massweight = .true.
176  if (.not.present(bathyt)) call mom_error(fatal, "int_density_dz_generic: "//&
177  "bathyT must be present if useMassWghtInterp is present and true.")
178  if (.not.present(dz_neglect)) call mom_error(fatal, "int_density_dz_generic: "//&
179  "dz_neglect must be present if useMassWghtInterp is present and true.")
180  endif ; endif
181 
182  do j=jsq,jeq+1 ; do i=isq,ieq+1
183  dz = z_t(i,j) - z_b(i,j)
184  do n=1,5
185  t5(n) = t(i,j) ; s5(n) = s(i,j)
186  p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
187  enddo
188  if (rho_scale /= 1.0) then
189  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
190  else
191  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks)
192  endif
193 
194  ! Use Boole's rule to estimate the pressure anomaly change.
195  rho_anom = c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))
196  dpa(i,j) = g_e*dz*rho_anom
197  ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
198  ! the pressure anomaly.
199  if (present(intz_dpa)) intz_dpa(i,j) = 0.5*g_e*dz**2 * &
200  (rho_anom - c1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
201  enddo ; enddo
202 
203  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
204  ! hWght is the distance measure by which the cell is violation of
205  ! hydrostatic consistency. For large hWght we bias the interpolation of
206  ! T & S along the top and bottom integrals, akin to thickness weighting.
207  hwght = 0.0
208  if (do_massweight) &
209  hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
210  if (hwght > 0.) then
211  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
212  hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
213  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
214  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
215  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
216  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
217  else
218  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
219  endif
220 
221  intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
222  do m=2,4
223  ! T, S, and z are interpolated in the horizontal. The z interpolation
224  ! is linear, but for T and S it may be thickness weighted.
225  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
226  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
227  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
228  t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
229  s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
230  p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i+1,j))
231  do n=2,5
232  t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) + gxrho*0.25*dz
233  enddo
234  if (rho_scale /= 1.0) then
235  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
236  else
237  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks)
238  endif
239 
240  ! Use Boole's rule to estimate the pressure anomaly change.
241  intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
242  enddo
243  ! Use Boole's rule to integrate the bottom pressure anomaly values in x.
244  intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
245  12.0*intz(3))
246  enddo ; enddo ; endif
247 
248  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
249  ! hWght is the distance measure by which the cell is violation of
250  ! hydrostatic consistency. For large hWght we bias the interpolation of
251  ! T & S along the top and bottom integrals, akin to thickness weighting.
252  hwght = 0.0
253  if (do_massweight) &
254  hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
255  if (hwght > 0.) then
256  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
257  hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
258  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
259  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
260  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
261  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
262  else
263  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
264  endif
265 
266  intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
267  do m=2,4
268  ! T, S, and z are interpolated in the horizontal. The z interpolation
269  ! is linear, but for T and S it may be thickness weighted.
270  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
271  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
272  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
273  t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
274  s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
275  p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i,j+1))
276  do n=2,5
277  t5(n) = t5(1) ; s5(n) = s5(1)
278  p5(n) = p5(n-1) + gxrho*0.25*dz
279  enddo
280  if (rho_scale /= 1.0) then
281  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
282  else
283  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks)
284  endif
285 
286  ! Use Boole's rule to estimate the pressure anomaly change.
287  intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
288  enddo
289  ! Use Boole's rule to integrate the values.
290  inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
291  12.0*intz(3))
292  enddo ; enddo ; endif
293 end subroutine int_density_dz_generic_pcm
294 
295 
296 !> Compute pressure gradient force integrals by quadrature for the case where
297 !! T and S are linear profiles.
298 subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
299  rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, dpa, &
300  intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp)
301  integer, intent(in) :: k !< Layer index to calculate integrals for
302  type(hor_index_type), intent(in) :: hi !< Ocean horizontal index structures for the input arrays
303  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
304  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
305  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
306  intent(in) :: t_t !< Potential temperature at the cell top [degC]
307  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
308  intent(in) :: t_b !< Potential temperature at the cell bottom [degC]
309  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
310  intent(in) :: s_t !< Salinity at the cell top [ppt]
311  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
312  intent(in) :: s_b !< Salinity at the cell bottom [ppt]
313  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &
314  intent(in) :: e !< Height of interfaces [Z ~> m]
315  real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is subtracted
316  !! out to reduce the magnitude of each of the integrals.
317  real, intent(in) :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used to calculate
318  !! the pressure (as p~=-z*rho_0*G_e) used in the equation of state.
319  real, intent(in) :: g_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
320  real, intent(in) :: dz_subroundoff !< A minuscule thickness change [Z ~> m]
321  real, dimension(SZI_(HI),SZJ_(HI)), &
322  intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]
323  type(eos_type), pointer :: eos !< Equation of state structure
324  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
325  real, dimension(SZI_(HI),SZJ_(HI)), &
326  intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]
327  real, dimension(SZI_(HI),SZJ_(HI)), &
328  optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of
329  !! the pressure anomaly relative to the anomaly at the
330  !! top of the layer [R L2 Z T-2 ~> Pa Z]
331  real, dimension(SZIB_(HI),SZJ_(HI)), &
332  optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
333  !! pressure anomaly at the top and bottom of the layer
334  !! divided by the x grid spacing [R L2 T-2 ~> Pa]
335  real, dimension(SZI_(HI),SZJB_(HI)), &
336  optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
337  !! pressure anomaly at the top and bottom of the layer
338  !! divided by the y grid spacing [R L2 T-2 ~> Pa]
339  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
340  !! interpolate T/S for top and bottom integrals.
341 
342 ! This subroutine calculates (by numerical quadrature) integrals of
343 ! pressure anomalies across layers, which are required for calculating the
344 ! finite-volume form pressure accelerations in a Boussinesq model. The one
345 ! potentially dodgy assumption here is that rho_0 is used both in the denominator
346 ! of the accelerations, and in the pressure used to calculated density (the
347 ! latter being -z*rho_0*G_e). These two uses could be separated if need be.
348 !
349 ! It is assumed that the salinity and temperature profiles are linear in the
350 ! vertical. The top and bottom values within each layer are provided and
351 ! a linear interpolation is used to compute intermediate values.
352 
353  ! Local variables
354  real :: t5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Temperatures along a line of subgrid locations [degC]
355  real :: s5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Salinities along a line of subgrid locations [ppt]
356  real :: t25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temperature variance along a line of subgrid locations [degC2]
357  real :: ts5((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS temp-salt covariance along a line of subgrid locations [degC ppt]
358  real :: s25((5*hi%iscb+1):(5*(hi%iecb+2))) ! SGS salinity variance along a line of subgrid locations [ppt2]
359  real :: p5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Pressures along a line of subgrid locations, never
360  ! rescaled from Pa [Pa]
361  real :: r5((5*hi%iscb+1):(5*(hi%iecb+2))) ! Densities anomalies along a line of subgrid
362  ! locations [R ~> kg m-3] or [kg m-3]
363  real :: t15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Temperatures at an array of subgrid locations [degC]
364  real :: s15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Salinities at an array of subgrid locations [ppt]
365  real :: t215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temperature variance along a line of subgrid locations [degC2]
366  real :: ts15((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS temp-salt covariance along a line of subgrid locations [degC ppt]
367  real :: s215((15*hi%iscb+1):(15*(hi%iecb+1))) ! SGS salinity variance along a line of subgrid locations [ppt2]
368  real :: p15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Pressures at an array of subgrid locations [Pa]
369  real :: r15((15*hi%iscb+1):(15*(hi%iecb+1))) ! Densities at an array of subgrid locations
370  ! [R ~> kg m-3] or [kg m-3]
371  real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]
372  real :: rho_anom ! A density anomaly [R ~> kg m-3] or [kg m-3]
373  real :: w_left, w_right ! Left and right weights [nondim]
374  real :: intz(5) ! The gravitational acceleration times the integrals of density
375  ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]
376  real, parameter :: c1_90 = 1.0/90.0 ! A rational constant [nondim]
377  real :: gxrho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
378  real :: i_rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]
379  real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]
380  real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]
381  real :: dz(hi%iscb:hi%iecb+1) ! Layer thicknesses at tracer points [Z ~> m]
382  real :: dz_x(5,hi%iscb:hi%iecb) ! Layer thicknesses along an x-line of subrid locations [Z ~> m]
383  real :: dz_y(5,hi%isc:hi%iec) ! Layer thicknesses along a y-line of subrid locations [Z ~> m]
384  real :: massweighttoggle ! A non-dimensional toggle factor (0 or 1) [nondim]
385  real :: ttl, tbl, ttr, tbr ! Temperatures at the velocity cell corners [degC]
386  real :: stl, sbl, str, sbr ! Salinities at the velocity cell corners [ppt]
387  real :: hwght ! A topographically limited thicknes weight [Z ~> m]
388  real :: hl, hr ! Thicknesses to the left and right [Z ~> m]
389  real :: idenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
390  logical :: use_stanley_eos ! True is SGS variance fields exist in tv.
391  logical :: use_vart, use_vars, use_covarts
392  integer :: isq, ieq, jsq, jeq, i, j, m, n
393  integer :: pos
394 
395  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
396 
397  rho_scale = us%kg_m3_to_R
398  gxrho = us%RL2_T2_to_Pa * g_e * rho_0
399  rho_ref_mks = rho_ref * us%R_to_kg_m3
400  i_rho = 1.0 / rho_0
401  massweighttoggle = 0.
402  if (present(usemasswghtinterp)) then
403  if (usemasswghtinterp) massweighttoggle = 1.
404  endif
405 
406  use_vart = associated(tv%varT)
407  use_covarts = associated(tv%covarTS)
408  use_vars = associated(tv%varS)
409  use_stanley_eos = use_vart .or. use_covarts .or. use_vars
410  t25(:) = 0.
411  ts5(:) = 0.
412  s25(:) = 0.
413  t215(:) = 0.
414  ts15(:) = 0.
415  s215(:) = 0.
416 
417  do n = 1, 5
418  wt_t(n) = 0.25 * real(5-n)
419  wt_b(n) = 1.0 - wt_t(n)
420  enddo
421 
422  ! 1. Compute vertical integrals
423  do j=jsq,jeq+1
424  do i = isq,ieq+1
425  dz(i) = e(i,j,k) - e(i,j,k+1)
426  do n=1,5
427  p5(i*5+n) = -gxrho*(e(i,j,k) - 0.25*real(n-1)*dz(i))
428  ! Salinity and temperature points are linearly interpolated
429  s5(i*5+n) = wt_t(n) * s_t(i,j,k) + wt_b(n) * s_b(i,j,k)
430  t5(i*5+n) = wt_t(n) * t_t(i,j,k) + wt_b(n) * t_b(i,j,k)
431  enddo
432  if (use_vart) t25(i*5+1:i*5+5) = tv%varT(i,j,k)
433  if (use_covarts) ts5(i*5+1:i*5+5) = tv%covarTS(i,j,k)
434  if (use_vars) s25(i*5+1:i*5+5) = tv%varS(i,j,k)
435  enddo
436  if (use_stanley_eos) then
437  if (rho_scale /= 1.0) then
438  call calculate_density(t5, s5, p5, t25, ts5, s25, r5, 1, (ieq-isq+2)*5, eos, &
439  rho_ref=rho_ref_mks, scale=rho_scale)
440  else
441  call calculate_density(t5, s5, p5, t25, ts5, s25, r5, 1, (ieq-isq+2)*5, eos, &
442  rho_ref=rho_ref_mks)
443  endif
444  else
445  if (rho_scale /= 1.0) then
446  call calculate_density(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho_ref=rho_ref_mks, &
447  scale=rho_scale)
448  else
449  call calculate_density(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho_ref=rho_ref_mks)
450  endif
451  endif
452 
453  do i=isq,ieq+1
454  ! Use Boole's rule to estimate the pressure anomaly change.
455  rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3))
456  dpa(i,j) = g_e*dz(i)*rho_anom
457  if (present(intz_dpa)) then
458  ! Use a Boole's-rule-like fifth-order accurate estimate of
459  ! the double integral of the pressure anomaly.
460  intz_dpa(i,j) = 0.5*g_e*dz(i)**2 * &
461  (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )
462  endif
463  enddo
464  enddo ! end loops on j
465 
466  ! 2. Compute horizontal integrals in the x direction
467  if (present(intx_dpa)) then ; do j=hi%jsc,hi%jec
468  do i=isq,ieq
469  ! Corner values of T and S
470  ! hWght is the distance measure by which the cell is violation of
471  ! hydrostatic consistency. For large hWght we bias the interpolation
472  ! of T,S along the top and bottom integrals, almost like thickness
473  ! weighting.
474  ! Note: To work in terrain following coordinates we could offset
475  ! this distance by the layer thickness to replicate other models.
476  hwght = massweighttoggle * &
477  max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k))
478  if (hwght > 0.) then
479  hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
480  hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz_subroundoff
481  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
482  idenom = 1./( hwght*(hr + hl) + hl*hr )
483  ttl = ( (hwght*hr)*t_t(i+1,j,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
484  ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i+1,j,k) ) * idenom
485  tbl = ( (hwght*hr)*t_b(i+1,j,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
486  tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i+1,j,k) ) * idenom
487  stl = ( (hwght*hr)*s_t(i+1,j,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
488  str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i+1,j,k) ) * idenom
489  sbl = ( (hwght*hr)*s_b(i+1,j,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
490  sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i+1,j,k) ) * idenom
491  else
492  ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i+1,j,k); tbr = t_b(i+1,j,k)
493  stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i+1,j,k); sbr = s_b(i+1,j,k)
494  endif
495 
496  do m=2,4
497  w_left = wt_t(m) ; w_right = wt_b(m)
498  dz_x(m,i) = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i+1,j,k) - e(i+1,j,k+1))
499 
500  ! Salinity and temperature points are linearly interpolated in
501  ! the horizontal. The subscript (1) refers to the top value in
502  ! the vertical profile while subscript (5) refers to the bottom
503  ! value in the vertical profile.
504  pos = i*15+(m-2)*5
505  t15(pos+1) = w_left*ttl + w_right*ttr
506  t15(pos+5) = w_left*tbl + w_right*tbr
507 
508  s15(pos+1) = w_left*stl + w_right*str
509  s15(pos+5) = w_left*sbl + w_right*sbr
510 
511  p15(pos+1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i+1,j,k))
512 
513  ! Pressure
514  do n=2,5
515  p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)
516  enddo
517 
518  ! Salinity and temperature (linear interpolation in the vertical)
519  do n=2,4
520  s15(pos+n) = wt_t(n) * s15(pos+1) + wt_b(n) * s15(pos+5)
521  t15(pos+n) = wt_t(n) * t15(pos+1) + wt_b(n) * t15(pos+5)
522  enddo
523  if (use_vart) t215(pos+1:pos+5) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i+1,j,k)
524  if (use_covarts) ts15(pos+1:pos+5) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i+1,j,k)
525  if (use_vars) s215(pos+1:pos+5) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i+1,j,k)
526  enddo
527  enddo
528 
529  if (use_stanley_eos) then
530  if (rho_scale /= 1.0) then
531  call calculate_density(t15, s15, p15, t215, ts15, s215, r15, 1, 15*(ieq-isq+1), eos, &
532  rho_ref=rho_ref_mks, scale=rho_scale)
533  else
534  call calculate_density(t15, s15, p15, t215, ts15, s215, r15, 1, 15*(ieq-isq+1), eos, &
535  rho_ref=rho_ref_mks)
536  endif
537  else
538  if (rho_scale /= 1.0) then
539  call calculate_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho_ref=rho_ref_mks, &
540  scale=rho_scale)
541  else
542  call calculate_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho_ref=rho_ref_mks)
543  endif
544  endif
545 
546  do i=isq,ieq
547  intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
548 
549  ! Use Boole's rule to estimate the pressure anomaly change.
550  do m = 2,4
551  pos = i*15+(m-2)*5
552  intz(m) = g_e*dz_x(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &
553  12.0*r15(pos+3)))
554  enddo
555  ! Use Boole's rule to integrate the bottom pressure anomaly values in x.
556  intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
557  12.0*intz(3))
558  enddo
559  enddo ; endif
560 
561  ! 3. Compute horizontal integrals in the y direction
562  if (present(inty_dpa)) then ; do j=jsq,jeq
563  do i=hi%isc,hi%iec
564  ! Corner values of T and S
565  ! hWght is the distance measure by which the cell is violation of
566  ! hydrostatic consistency. For large hWght we bias the interpolation
567  ! of T,S along the top and bottom integrals, almost like thickness
568  ! weighting.
569  ! Note: To work in terrain following coordinates we could offset
570  ! this distance by the layer thickness to replicate other models.
571  hwght = massweighttoggle * &
572  max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k))
573  if (hwght > 0.) then
574  hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
575  hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz_subroundoff
576  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
577  idenom = 1./( hwght*(hr + hl) + hl*hr )
578  ttl = ( (hwght*hr)*t_t(i,j+1,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
579  ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i,j+1,k) ) * idenom
580  tbl = ( (hwght*hr)*t_b(i,j+1,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
581  tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i,j+1,k) ) * idenom
582  stl = ( (hwght*hr)*s_t(i,j+1,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
583  str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i,j+1,k) ) * idenom
584  sbl = ( (hwght*hr)*s_b(i,j+1,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
585  sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i,j+1,k) ) * idenom
586  else
587  ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i,j+1,k); tbr = t_b(i,j+1,k)
588  stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i,j+1,k); sbr = s_b(i,j+1,k)
589  endif
590 
591  do m=2,4
592  w_left = wt_t(m) ; w_right = wt_b(m)
593  dz_y(m,i) = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i,j+1,k) - e(i,j+1,k+1))
594 
595  ! Salinity and temperature points are linearly interpolated in
596  ! the horizontal. The subscript (1) refers to the top value in
597  ! the vertical profile while subscript (5) refers to the bottom
598  ! value in the vertical profile.
599  pos = i*15+(m-2)*5
600  t15(pos+1) = w_left*ttl + w_right*ttr
601  t15(pos+5) = w_left*tbl + w_right*tbr
602 
603  s15(pos+1) = w_left*stl + w_right*str
604  s15(pos+5) = w_left*sbl + w_right*sbr
605 
606  p15(pos+1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i,j+1,k))
607 
608  ! Pressure
609  do n=2,5
610  p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i)
611  enddo
612 
613  ! Salinity and temperature (linear interpolation in the vertical)
614  do n=2,4
615  s15(pos+n) = wt_t(n) * s15(pos+1) + wt_b(n) * s15(pos+5)
616  t15(pos+n) = wt_t(n) * t15(pos+1) + wt_b(n) * t15(pos+5)
617  enddo
618  if (use_vart) t215(pos+1:pos+5) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i,j+1,k)
619  if (use_covarts) ts15(pos+1:pos+5) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i,j+1,k)
620  if (use_vars) s215(pos+1:pos+5) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i,j+1,k)
621  enddo
622  enddo
623 
624  if (use_stanley_eos) then
625  if (rho_scale /= 1.0) then
626  call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
627  t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), &
628  r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, &
629  rho_ref=rho_ref_mks, scale=rho_scale)
630  else
631  call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
632  t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), &
633  r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, rho_ref=rho_ref_mks)
634  endif
635  else
636  if (rho_scale /= 1.0) then
637  call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
638  r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, &
639  rho_ref=rho_ref_mks, scale=rho_scale)
640  else
641  call calculate_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), &
642  r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, rho_ref=rho_ref_mks)
643  endif
644  endif
645 
646  do i=hi%isc,hi%iec
647  intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
648 
649  ! Use Boole's rule to estimate the pressure anomaly change.
650  do m = 2,4
651  pos = i*15+(m-2)*5
652  intz(m) = g_e*dz_y(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &
653  32.0*(r15(pos+2)+r15(pos+4)) + &
654  12.0*r15(pos+3)))
655  enddo
656  ! Use Boole's rule to integrate the values.
657  inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
658  12.0*intz(3))
659  enddo
660  enddo ; endif
661 
662 end subroutine int_density_dz_generic_plm
663 
664 
665 !> Compute pressure gradient force integrals for layer "k" and the case where T and S
666 !! are parabolic profiles
667 subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
668  rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, &
669  dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp)
670  integer, intent(in) :: k !< Layer index to calculate integrals for
671  type(hor_index_type), intent(in) :: hi !< Ocean horizontal index structures for the input arrays
672  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
673  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
674  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
675  intent(in) :: t_t !< Potential temperature at the cell top [degC]
676  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
677  intent(in) :: t_b !< Potential temperature at the cell bottom [degC]
678  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
679  intent(in) :: s_t !< Salinity at the cell top [ppt]
680  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
681  intent(in) :: s_b !< Salinity at the cell bottom [ppt]
682  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &
683  intent(in) :: e !< Height of interfaces [Z ~> m]
684  real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is
685  !! subtracted out to reduce the magnitude of each of the integrals.
686  real, intent(in) :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used to calculate
687  !! the pressure (as p~=-z*rho_0*G_e) used in the equation of state.
688  real, intent(in) :: g_e !< The Earth's gravitational acceleration [m s-2]
689  real, intent(in) :: dz_subroundoff !< A minuscule thickness change [Z ~> m]
690  real, dimension(SZI_(HI),SZJ_(HI)), &
691  intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]
692  type(eos_type), pointer :: eos !< Equation of state structure
693  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
694  real, dimension(SZI_(HI),SZJ_(HI)), &
695  intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]
696  real, dimension(SZI_(HI),SZJ_(HI)), &
697  optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of
698  !! the pressure anomaly relative to the anomaly at the
699  !! top of the layer [R L2 Z T-2 ~> Pa m]
700  real, dimension(SZIB_(HI),SZJ_(HI)), &
701  optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
702  !! pressure anomaly at the top and bottom of the layer
703  !! divided by the x grid spacing [R L2 T-2 ~> Pa]
704  real, dimension(SZI_(HI),SZJB_(HI)), &
705  optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
706  !! pressure anomaly at the top and bottom of the layer
707  !! divided by the y grid spacing [R L2 T-2 ~> Pa]
708  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
709  !! interpolate T/S for top and bottom integrals.
710 
711 ! This subroutine calculates (by numerical quadrature) integrals of
712 ! pressure anomalies across layers, which are required for calculating the
713 ! finite-volume form pressure accelerations in a Boussinesq model. The one
714 ! potentially dodgy assumption here is that rho_0 is used both in the denominator
715 ! of the accelerations, and in the pressure used to calculated density (the
716 ! latter being -z*rho_0*G_e). These two uses could be separated if need be.
717 !
718 ! It is assumed that the salinity and temperature profiles are parabolic in the
719 ! vertical. The top and bottom values within each layer are provided and
720 ! a parabolic interpolation is used to compute intermediate values.
721 
722  ! Local variables
723  real :: t5(5) ! Temperatures along a line of subgrid locations [degC]
724  real :: s5(5) ! Salinities along a line of subgrid locations [ppt]
725  real :: t25(5) ! SGS temperature variance along a line of subgrid locations [degC2]
726  real :: ts5(5) ! SGS temperature-salinity covariance along a line of subgrid locations [degC ppt]
727  real :: s25(5) ! SGS salinity variance along a line of subgrid locations [ppt2]
728  real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa]
729  real :: r5(5) ! Density anomalies from rho_ref at quadrature points [R ~> kg m-3] or [kg m-3]
730  real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]
731  real :: rho_anom ! The integrated density anomaly [R ~> kg m-3] or [kg m-3]
732  real :: w_left, w_right ! Left and right weights [nondim]
733  real :: intz(5) ! The gravitational acceleration times the integrals of density
734  ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]
735  real, parameter :: c1_90 = 1.0/90.0 ! Rational constants.
736  real :: gxrho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
737  real :: i_rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]
738  real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]
739  real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]
740  real :: dz ! Layer thicknesses at tracer points [Z ~> m]
741  real :: massweighttoggle ! A non-dimensional toggle factor (0 or 1) [nondim]
742  real :: ttl, tbl, tml, ttr, tbr, tmr ! Temperatures at the velocity cell corners [degC]
743  real :: stl, sbl, sml, str, sbr, smr ! Salinities at the velocity cell corners [ppt]
744  real :: s6 ! PPM curvature coefficient for S [ppt]
745  real :: t6 ! PPM curvature coefficient for T [degC]
746  real :: t_top, t_mn, t_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T
747  real :: s_top, s_mn, s_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S
748  real :: hwght ! A topographically limited thicknes weight [Z ~> m]
749  real :: hl, hr ! Thicknesses to the left and right [Z ~> m]
750  real :: idenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
751  integer :: isq, ieq, jsq, jeq, i, j, m, n
752  logical :: use_ppm ! If false, assume zero curvature in reconstruction, i.e. PLM
753  logical :: use_stanley_eos ! True is SGS variance fields exist in tv.
754  logical :: use_vart, use_vars, use_covarts
755 
756  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
757 
758  rho_scale = us%kg_m3_to_R
759  gxrho = us%RL2_T2_to_Pa * g_e * rho_0
760  rho_ref_mks = rho_ref * us%R_to_kg_m3
761  i_rho = 1.0 / rho_0
762  massweighttoggle = 0.
763  if (present(usemasswghtinterp)) then
764  if (usemasswghtinterp) massweighttoggle = 1.
765  endif
766 
767  ! In event PPM calculation is bypassed with use_PPM=False
768  s6 = 0.
769  t6 = 0.
770  use_ppm = .true. ! This is a place-holder to allow later re-use of this function
771 
772  use_vart = associated(tv%varT)
773  use_covarts = associated(tv%covarTS)
774  use_vars = associated(tv%varS)
775  use_stanley_eos = use_vart .or. use_covarts .or. use_vars
776  t25(:) = 0.
777  ts5(:) = 0.
778  s25(:) = 0.
779 
780  do n = 1, 5
781  wt_t(n) = 0.25 * real(5-n)
782  wt_b(n) = 1.0 - wt_t(n)
783  enddo
784 
785  ! 1. Compute vertical integrals
786  do j=jsq,jeq+1 ; do i=isq,ieq+1
787  if (use_ppm) then
788  ! Curvature coefficient of the parabolas
789  s6 = 3.0 * ( 2.0*tv%S(i,j,k) - ( s_t(i,j,k) + s_b(i,j,k) ) )
790  t6 = 3.0 * ( 2.0*tv%T(i,j,k) - ( t_t(i,j,k) + t_b(i,j,k) ) )
791  endif
792  dz = e(i,j,k) - e(i,j,k+1)
793  do n=1,5
794  p5(n) = -gxrho*(e(i,j,k) - 0.25*real(n-1)*dz)
795  ! Salinity and temperature points are reconstructed with PPM
796  s5(n) = wt_t(n) * s_t(i,j,k) + wt_b(n) * ( s_b(i,j,k) + s6 * wt_t(n) )
797  t5(n) = wt_t(n) * t_t(i,j,k) + wt_b(n) * ( t_b(i,j,k) + t6 * wt_t(n) )
798  enddo
799  if (use_stanley_eos) then
800  if (use_vart) t25(:) = tv%varT(i,j,k)
801  if (use_covarts) ts5(:) = tv%covarTS(i,j,k)
802  if (use_vars) s25(:) = tv%varS(i,j,k)
803  call calculate_density(t5, s5, p5, t25, ts5, s25, r5, &
804  1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
805  else
806  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
807  endif
808 
809  ! Use Boole's rule to estimate the pressure anomaly change.
810  rho_anom = c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))
811  dpa(i,j) = g_e*dz*rho_anom
812  if (present(intz_dpa)) then
813  ! Use a Boole's-rule-like fifth-order accurate estimate of
814  ! the double integral of the pressure anomaly.
815  intz_dpa(i,j) = 0.5*g_e*dz**2 * &
816  (rho_anom - c1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
817  endif
818  enddo ; enddo ! end loops on j and i
819 
820  ! 2. Compute horizontal integrals in the x direction
821  if (present(intx_dpa)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
822  ! Corner values of T and S
823  ! hWght is the distance measure by which the cell is violation of
824  ! hydrostatic consistency. For large hWght we bias the interpolation
825  ! of T,S along the top and bottom integrals, almost like thickness
826  ! weighting.
827  ! Note: To work in terrain following coordinates we could offset
828  ! this distance by the layer thickness to replicate other models.
829  hwght = massweighttoggle * &
830  max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k))
831  if (hwght > 0.) then
832  hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
833  hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz_subroundoff
834  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
835  idenom = 1./( hwght*(hr + hl) + hl*hr )
836  ttl = ( (hwght*hr)*t_t(i+1,j,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
837  tbl = ( (hwght*hr)*t_b(i+1,j,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
838  tml = ( (hwght*hr)*tv%T(i+1,j,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom
839  ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i+1,j,k) ) * idenom
840  tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i+1,j,k) ) * idenom
841  tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i+1,j,k) ) * idenom
842  stl = ( (hwght*hr)*s_t(i+1,j,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
843  sbl = ( (hwght*hr)*s_b(i+1,j,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
844  sml = ( (hwght*hr)*tv%S(i+1,j,k) + (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom
845  str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i+1,j,k) ) * idenom
846  sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i+1,j,k) ) * idenom
847  smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i+1,j,k) ) * idenom
848  else
849  ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i+1,j,k); tbr = t_b(i+1,j,k)
850  tml = tv%T(i,j,k); tmr = tv%T(i+1,j,k)
851  stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i+1,j,k); sbr = s_b(i+1,j,k)
852  sml = tv%S(i,j,k); smr = tv%S(i+1,j,k)
853  endif
854 
855  do m=2,4
856  w_left = wt_t(m) ; w_right = wt_b(m)
857 
858  ! Salinity and temperature points are linearly interpolated in
859  ! the horizontal. The subscript (1) refers to the top value in
860  ! the vertical profile while subscript (5) refers to the bottom
861  ! value in the vertical profile.
862  t_top = w_left*ttl + w_right*ttr
863  t_mn = w_left*tml + w_right*tmr
864  t_bot = w_left*tbl + w_right*tbr
865 
866  s_top = w_left*stl + w_right*str
867  s_mn = w_left*sml + w_right*smr
868  s_bot = w_left*sbl + w_right*sbr
869 
870  ! Pressure
871  dz = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i+1,j,k) - e(i+1,j,k+1))
872  p5(1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i+1,j,k))
873  do n=2,5
874  p5(n) = p5(n-1) + gxrho*0.25*dz
875  enddo
876 
877  ! Parabolic reconstructions in the vertical for T and S
878  if (use_ppm) then
879  ! Coefficients of the parabolas
880  s6 = 3.0 * ( 2.0*s_mn - ( s_top + s_bot ) )
881  t6 = 3.0 * ( 2.0*t_mn - ( t_top + t_bot ) )
882  endif
883  do n=1,5
884  s5(n) = wt_t(n) * s_top + wt_b(n) * ( s_bot + s6 * wt_t(n) )
885  t5(n) = wt_t(n) * t_top + wt_b(n) * ( t_bot + t6 * wt_t(n) )
886  enddo
887 
888  if (use_stanley_eos) then
889  if (use_vart) t25(:) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i+1,j,k)
890  if (use_covarts) ts5(:) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i+1,j,k)
891  if (use_vars) s25(:) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i+1,j,k)
892  call calculate_density(t5, s5, p5, t25, ts5, s25, r5, &
893  1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
894  else
895  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
896  endif
897 
898  ! Use Boole's rule to estimate the pressure anomaly change.
899  intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) )
900  enddo ! m
901  intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
902 
903  ! Use Boole's rule to integrate the bottom pressure anomaly values in x.
904  intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))
905 
906  enddo ; enddo ; endif
907 
908  ! 3. Compute horizontal integrals in the y direction
909  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
910  ! Corner values of T and S
911  ! hWght is the distance measure by which the cell is violation of
912  ! hydrostatic consistency. For large hWght we bias the interpolation
913  ! of T,S along the top and bottom integrals, almost like thickness
914  ! weighting.
915  ! Note: To work in terrain following coordinates we could offset
916  ! this distance by the layer thickness to replicate other models.
917  hwght = massweighttoggle * &
918  max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k))
919  if (hwght > 0.) then
920  hl = (e(i,j,k) - e(i,j,k+1)) + dz_subroundoff
921  hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz_subroundoff
922  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
923  idenom = 1./( hwght*(hr + hl) + hl*hr )
924  ttl = ( (hwght*hr)*t_t(i,j+1,k) + (hwght*hl + hr*hl)*t_t(i,j,k) ) * idenom
925  tbl = ( (hwght*hr)*t_b(i,j+1,k) + (hwght*hl + hr*hl)*t_b(i,j,k) ) * idenom
926  tml = ( (hwght*hr)*tv%T(i,j+1,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom
927  ttr = ( (hwght*hl)*t_t(i,j,k) + (hwght*hr + hr*hl)*t_t(i,j+1,k) ) * idenom
928  tbr = ( (hwght*hl)*t_b(i,j,k) + (hwght*hr + hr*hl)*t_b(i,j+1,k) ) * idenom
929  tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i,j+1,k) ) * idenom
930  stl = ( (hwght*hr)*s_t(i,j+1,k) + (hwght*hl + hr*hl)*s_t(i,j,k) ) * idenom
931  sbl = ( (hwght*hr)*s_b(i,j+1,k) + (hwght*hl + hr*hl)*s_b(i,j,k) ) * idenom
932  sml = ( (hwght*hr)*tv%S(i,j+1,k)+ (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom
933  str = ( (hwght*hl)*s_t(i,j,k) + (hwght*hr + hr*hl)*s_t(i,j+1,k) ) * idenom
934  sbr = ( (hwght*hl)*s_b(i,j,k) + (hwght*hr + hr*hl)*s_b(i,j+1,k) ) * idenom
935  smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i,j+1,k) ) * idenom
936  else
937  ttl = t_t(i,j,k); tbl = t_b(i,j,k); ttr = t_t(i,j+1,k); tbr = t_b(i,j+1,k)
938  tml = tv%T(i,j,k); tmr = tv%T(i,j+1,k)
939  stl = s_t(i,j,k); sbl = s_b(i,j,k); str = s_t(i,j+1,k); sbr = s_b(i,j+1,k)
940  sml = tv%S(i,j,k); smr = tv%S(i,j+1,k)
941  endif
942 
943  do m=2,4
944  w_left = wt_t(m) ; w_right = wt_b(m)
945 
946  ! Salinity and temperature points are linearly interpolated in
947  ! the horizontal. The subscript (1) refers to the top value in
948  ! the vertical profile while subscript (5) refers to the bottom
949  ! value in the vertical profile.
950  t_top = w_left*ttl + w_right*ttr
951  t_mn = w_left*tml + w_right*tmr
952  t_bot = w_left*tbl + w_right*tbr
953 
954  s_top = w_left*stl + w_right*str
955  s_mn = w_left*sml + w_right*smr
956  s_bot = w_left*sbl + w_right*sbr
957 
958  ! Pressure
959  dz = w_left*(e(i,j,k) - e(i,j,k+1)) + w_right*(e(i,j+1,k) - e(i,j+1,k+1))
960  p5(1) = -gxrho*(w_left*e(i,j,k) + w_right*e(i,j+1,k))
961  do n=2,5
962  p5(n) = p5(n-1) + gxrho*0.25*dz
963  enddo
964 
965  ! Parabolic reconstructions in the vertical for T and S
966  if (use_ppm) then
967  ! Coefficients of the parabolas
968  s6 = 3.0 * ( 2.0*s_mn - ( s_top + s_bot ) )
969  t6 = 3.0 * ( 2.0*t_mn - ( t_top + t_bot ) )
970  endif
971  do n=1,5
972  s5(n) = wt_t(n) * s_top + wt_b(n) * ( s_bot + s6 * wt_t(n) )
973  t5(n) = wt_t(n) * t_top + wt_b(n) * ( t_bot + t6 * wt_t(n) )
974  enddo
975 
976  if (use_stanley_eos) then
977  if (use_vart) t25(:) = w_left*tv%varT(i,j,k) + w_right*tv%varT(i,j+1,k)
978  if (use_covarts) ts5(:) = w_left*tv%covarTS(i,j,k) + w_right*tv%covarTS(i,j+1,k)
979  if (use_vars) s25(:) = w_left*tv%varS(i,j,k) + w_right*tv%varS(i,j+1,k)
980  call calculate_density(t5, s5, p5, t25, ts5, s25, r5, &
981  1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
982  else
983  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref=rho_ref_mks, scale=rho_scale)
984  endif
985 
986  ! Use Boole's rule to estimate the pressure anomaly change.
987  intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) )
988  enddo ! m
989  intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
990 
991  ! Use Boole's rule to integrate the bottom pressure anomaly values in y.
992  inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))
993 
994  enddo ; enddo ; endif
995 
996 end subroutine int_density_dz_generic_ppm
997 
998 !> Calls the appropriate subroutine to calculate analytical and nearly-analytical
999 !! integrals in pressure across layers of geopotential anomalies, which are
1000 !! required for calculating the finite-volume form pressure accelerations in a
1001 !! non-Boussinesq model. There are essentially no free assumptions, apart from the
1002 !! use of Boole's rule to do the horizontal integrals, and from a truncation in the
1003 !! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
1004 subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, &
1005  dza, intp_dza, intx_dza, inty_dza, halo_size, &
1006  bathyP, dP_tiny, useMassWghtInterp)
1007  type(hor_index_type), intent(in) :: hi !< The horizontal index structure
1008  real, dimension(SZI_(HI),SZJ_(HI)), &
1009  intent(in) :: t !< Potential temperature referenced to the surface [degC]
1010  real, dimension(SZI_(HI),SZJ_(HI)), &
1011  intent(in) :: s !< Salinity [ppt]
1012  real, dimension(SZI_(HI),SZJ_(HI)), &
1013  intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]
1014  real, dimension(SZI_(HI),SZJ_(HI)), &
1015  intent(in) :: p_b !< Pressure at the bottom of the layer [R L2 T-2 ~> Pa] or [Pa]
1016  real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
1017  !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
1018  !! The calculation is mathematically identical with different values of
1019  !! alpha_ref, but this reduces the effects of roundoff.
1020  type(eos_type), pointer :: eos !< Equation of state structure
1021  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1022  real, dimension(SZI_(HI),SZJ_(HI)), &
1023  intent(inout) :: dza !< The change in the geopotential anomaly across
1024  !! the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]
1025  real, dimension(SZI_(HI),SZJ_(HI)), &
1026  optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the
1027  !! geopotential anomaly relative to the anomaly at the bottom of the
1028  !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]
1029  real, dimension(SZIB_(HI),SZJ_(HI)), &
1030  optional, intent(inout) :: intx_dza !< The integral in x of the difference between the
1031  !! geopotential anomaly at the top and bottom of the layer divided by
1032  !! the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
1033  real, dimension(SZI_(HI),SZJB_(HI)), &
1034  optional, intent(inout) :: inty_dza !< The integral in y of the difference between the
1035  !! geopotential anomaly at the top and bottom of the layer divided by
1036  !! the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
1037  integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
1038  real, dimension(SZI_(HI),SZJ_(HI)), &
1039  optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
1040  real, optional, intent(in) :: dp_tiny !< A minuscule pressure change with
1041  !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
1042  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
1043  !! to interpolate T/S for top and bottom integrals.
1044 
1045  if (eos_quadrature(eos)) then
1046  call int_spec_vol_dp_generic_pcm(t, s, p_t, p_b, alpha_ref, hi, eos, us, &
1047  dza, intp_dza, intx_dza, inty_dza, halo_size, &
1048  bathyp, dp_tiny, usemasswghtinterp)
1049  else
1050  call analytic_int_specific_vol_dp(t, s, p_t, p_b, alpha_ref, hi, eos, &
1051  dza, intp_dza, intx_dza, inty_dza, halo_size, &
1052  bathyp, dp_tiny, usemasswghtinterp)
1053  endif
1054 
1055 end subroutine int_specific_vol_dp
1056 
1057 
1058 !> This subroutine calculates integrals of specific volume anomalies in
1059 !! pressure across layers, which are required for calculating the finite-volume
1060 !! form pressure accelerations in a non-Boussinesq model. There are essentially
1061 !! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.
1062 subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, dza, &
1063  intp_dza, intx_dza, inty_dza, halo_size, &
1064  bathyP, dP_neglect, useMassWghtInterp)
1065  type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
1066  real, dimension(SZI_(HI),SZJ_(HI)), &
1067  intent(in) :: t !< Potential temperature of the layer [degC]
1068  real, dimension(SZI_(HI),SZJ_(HI)), &
1069  intent(in) :: s !< Salinity of the layer [ppt]
1070  real, dimension(SZI_(HI),SZJ_(HI)), &
1071  intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]
1072  real, dimension(SZI_(HI),SZJ_(HI)), &
1073  intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]
1074  real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
1075  !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
1076  !! The calculation is mathematically identical with different values of
1077  !! alpha_ref, but alpha_ref alters the effects of roundoff, and
1078  !! answers do change.
1079  type(eos_type), pointer :: eos !< Equation of state structure
1080  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1081  real, dimension(SZI_(HI),SZJ_(HI)), &
1082  intent(inout) :: dza !< The change in the geopotential anomaly
1083  !! across the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]
1084  real, dimension(SZI_(HI),SZJ_(HI)), &
1085  optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
1086  !! the geopotential anomaly relative to the anomaly at the bottom of the
1087  !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]
1088  real, dimension(SZIB_(HI),SZJ_(HI)), &
1089  optional, intent(inout) :: intx_dza !< The integral in x of the difference between
1090  !! the geopotential anomaly at the top and bottom of the layer divided
1091  !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
1092  real, dimension(SZI_(HI),SZJB_(HI)), &
1093  optional, intent(inout) :: inty_dza !< The integral in y of the difference between
1094  !! the geopotential anomaly at the top and bottom of the layer divided
1095  !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
1096  integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
1097  real, dimension(SZI_(HI),SZJ_(HI)), &
1098  optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
1099  real, optional, intent(in) :: dp_neglect !< A minuscule pressure change with
1100  !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
1101  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
1102  !! to interpolate T/S for top and bottom integrals.
1103 
1104 ! This subroutine calculates analytical and nearly-analytical integrals in
1105 ! pressure across layers of geopotential anomalies, which are required for
1106 ! calculating the finite-volume form pressure accelerations in a non-Boussinesq
1107 ! model. There are essentially no free assumptions, apart from the use of
1108 ! Boole's rule to do the horizontal integrals, and from a truncation in the
1109 ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
1110 
1111  ! Local variables
1112  real :: t5(5) ! Temperatures at five quadrature points [degC]
1113  real :: s5(5) ! Salinities at five quadrature points [ppt]
1114  real :: p5(5) ! Pressures at five quadrature points, scaled back to Pa if necessary [Pa]
1115  real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]
1116  real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]
1117  real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]
1118  real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
1119  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
1120  real :: alpha_ref_mks ! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1]
1121  real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
1122  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]
1123  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]
1124  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]
1125  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]
1126  real :: intp(5) ! The integrals of specific volume with pressure at the
1127  ! 5 sub-column locations [L2 T-2 ~> m2 s-2]
1128  real :: rl2_t2_to_pa ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1]
1129  real :: sv_scale ! A multiplicative factor by which to scale specific
1130  ! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
1131  logical :: do_massweight ! Indicates whether to do mass weighting.
1132  real, parameter :: c1_90 = 1.0/90.0 ! A rational constant.
1133  integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, n, halo
1134 
1135  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
1136  halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
1137  ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
1138  if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh); endif
1139  if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh); endif
1140 
1141  sv_scale = us%R_to_kg_m3
1142  rl2_t2_to_pa = us%RL2_T2_to_Pa
1143  alpha_ref_mks = alpha_ref * us%kg_m3_to_R
1144 
1145  do_massweight = .false.
1146  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
1147  do_massweight = .true.
1148  if (.not.present(bathyp)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
1149  "bathyP must be present if useMassWghtInterp is present and true.")
1150  if (.not.present(dp_neglect)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
1151  "dP_neglect must be present if useMassWghtInterp is present and true.")
1152  endif ; endif
1153 
1154  do j=jsh,jeh ; do i=ish,ieh
1155  dp = p_b(i,j) - p_t(i,j)
1156  do n=1,5
1157  t5(n) = t(i,j) ; s5(n) = s(i,j)
1158  p5(n) = rl2_t2_to_pa * (p_b(i,j) - 0.25*real(n-1)*dp)
1159  enddo
1160 
1161  if (sv_scale /= 1.0) then
1162  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
1163  else
1164  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
1165  endif
1166 
1167  ! Use Boole's rule to estimate the interface height anomaly change.
1168  alpha_anom = c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3))
1169  dza(i,j) = dp*alpha_anom
1170  ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
1171  ! the interface height anomaly.
1172  if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
1173  (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
1174  enddo ; enddo
1175 
1176  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
1177  ! hWght is the distance measure by which the cell is violation of
1178  ! hydrostatic consistency. For large hWght we bias the interpolation of
1179  ! T & S along the top and bottom integrals, akin to thickness weighting.
1180  hwght = 0.0
1181  if (do_massweight) &
1182  hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
1183  if (hwght > 0.) then
1184  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
1185  hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
1186  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1187  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1188  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1189  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1190  else
1191  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1192  endif
1193 
1194  intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
1195  do m=2,4
1196  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1197  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1198 
1199  ! T, S, and p are interpolated in the horizontal. The p interpolation
1200  ! is linear, but for T and S it may be thickness weighted.
1201  p5(1) = rl2_t2_to_pa * (wt_l*p_b(i,j) + wt_r*p_b(i+1,j))
1202  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
1203  t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
1204  s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
1205 
1206  do n=2,5
1207  t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - rl2_t2_to_pa * 0.25*dp
1208  enddo
1209  if (sv_scale /= 1.0) then
1210  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
1211  else
1212  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
1213  endif
1214 
1215  ! Use Boole's rule to estimate the interface height anomaly change.
1216  intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
1217  12.0*a5(3)))
1218  enddo
1219  ! Use Boole's rule to integrate the interface height anomaly values in x.
1220  intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
1221  12.0*intp(3))
1222  enddo ; enddo ; endif
1223 
1224  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
1225  ! hWght is the distance measure by which the cell is violation of
1226  ! hydrostatic consistency. For large hWght we bias the interpolation of
1227  ! T & S along the top and bottom integrals, akin to thickness weighting.
1228  hwght = 0.0
1229  if (do_massweight) &
1230  hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
1231  if (hwght > 0.) then
1232  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
1233  hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
1234  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1235  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1236  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1237  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1238  else
1239  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1240  endif
1241 
1242  intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
1243  do m=2,4
1244  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1245  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1246 
1247  ! T, S, and p are interpolated in the horizontal. The p interpolation
1248  ! is linear, but for T and S it may be thickness weighted.
1249  p5(1) = rl2_t2_to_pa * (wt_l*p_b(i,j) + wt_r*p_b(i,j+1))
1250  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
1251  t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
1252  s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
1253  do n=2,5
1254  t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = rl2_t2_to_pa * (p5(n-1) - 0.25*dp)
1255  enddo
1256  if (sv_scale /= 1.0) then
1257  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
1258  else
1259  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
1260  endif
1261 
1262  ! Use Boole's rule to estimate the interface height anomaly change.
1263  intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
1264  12.0*a5(3)))
1265  enddo
1266  ! Use Boole's rule to integrate the interface height anomaly values in y.
1267  inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
1268  12.0*intp(3))
1269  enddo ; enddo ; endif
1270 
1271 end subroutine int_spec_vol_dp_generic_pcm
1272 
1273 !> This subroutine calculates integrals of specific volume anomalies in
1274 !! pressure across layers, which are required for calculating the finite-volume
1275 !! form pressure accelerations in a non-Boussinesq model. There are essentially
1276 !! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.
1277 subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, &
1278  dP_neglect, bathyP, HI, EOS, US, dza, &
1279  intp_dza, intx_dza, inty_dza, useMassWghtInterp)
1280  type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
1281  real, dimension(SZI_(HI),SZJ_(HI)), &
1282  intent(in) :: t_t !< Potential temperature at the top of the layer [degC]
1283  real, dimension(SZI_(HI),SZJ_(HI)), &
1284  intent(in) :: t_b !< Potential temperature at the bottom of the layer [degC]
1285  real, dimension(SZI_(HI),SZJ_(HI)), &
1286  intent(in) :: s_t !< Salinity at the top the layer [ppt]
1287  real, dimension(SZI_(HI),SZJ_(HI)), &
1288  intent(in) :: s_b !< Salinity at the bottom the layer [ppt]
1289  real, dimension(SZI_(HI),SZJ_(HI)), &
1290  intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]
1291  real, dimension(SZI_(HI),SZJ_(HI)), &
1292  intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]
1293  real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
1294  !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
1295  !! The calculation is mathematically identical with different values of
1296  !! alpha_ref, but alpha_ref alters the effects of roundoff, and
1297  !! answers do change.
1298  real, intent(in) :: dp_neglect !<!< A miniscule pressure change with
1299  !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
1300  real, dimension(SZI_(HI),SZJ_(HI)), &
1301  intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
1302  type(eos_type), pointer :: eos !< Equation of state structure
1303  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1304  real, dimension(SZI_(HI),SZJ_(HI)), &
1305  intent(inout) :: dza !< The change in the geopotential anomaly
1306  !! across the layer [L2 T-2 ~> m2 s-2]
1307  real, dimension(SZI_(HI),SZJ_(HI)), &
1308  optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
1309  !! the geopotential anomaly relative to the anomaly at the bottom of the
1310  !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]
1311  real, dimension(SZIB_(HI),SZJ_(HI)), &
1312  optional, intent(inout) :: intx_dza !< The integral in x of the difference between
1313  !! the geopotential anomaly at the top and bottom of the layer divided
1314  !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
1315  real, dimension(SZI_(HI),SZJB_(HI)), &
1316  optional, intent(inout) :: inty_dza !< The integral in y of the difference between
1317  !! the geopotential anomaly at the top and bottom of the layer divided
1318  !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]
1319  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
1320  !! to interpolate T/S for top and bottom integrals.
1321 
1322 ! This subroutine calculates analytical and nearly-analytical integrals in
1323 ! pressure across layers of geopotential anomalies, which are required for
1324 ! calculating the finite-volume form pressure accelerations in a non-Boussinesq
1325 ! model. There are essentially no free assumptions, apart from the use of
1326 ! Boole's rule to do the horizontal integrals, and from a truncation in the
1327 ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
1328 
1329  real :: t5(5) ! Temperatures at five quadrature points [degC]
1330  real :: s5(5) ! Salinities at five quadrature points [ppt]
1331  real :: p5(5) ! Pressures at five quadrature points, scaled back to Pa as necessary [Pa]
1332  real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]
1333  real :: t15(15) ! Temperatures at fifteen interior quadrature points [degC]
1334  real :: s15(15) ! Salinities at fifteen interior quadrature points [ppt]
1335  real :: p15(15) ! Pressures at fifteen quadrature points, scaled back to Pa as necessary [Pa]
1336  real :: a15(15) ! Specific volumes at fifteen quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]
1337  real :: wt_t(5), wt_b(5) ! Weights of top and bottom values at quadrature points [nondim]
1338  real :: t_top, t_bot, s_top, s_bot, p_top, p_bot
1339 
1340  real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1]
1341  real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]
1342  real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa]
1343  real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
1344  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
1345  real :: alpha_ref_mks ! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1]
1346  real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
1347  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim]
1348  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim]
1349  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim]
1350  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim]
1351  real :: intp(5) ! The integrals of specific volume with pressure at the
1352  ! 5 sub-column locations [L2 T-2 ~> m2 s-2]
1353  real :: rl2_t2_to_pa ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1]
1354  real :: sv_scale ! A multiplicative factor by which to scale specific
1355  ! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
1356  real, parameter :: c1_90 = 1.0/90.0 ! A rational constant.
1357  logical :: do_massweight ! Indicates whether to do mass weighting.
1358  integer :: isq, ieq, jsq, jeq, i, j, m, n, pos
1359 
1360  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
1361 
1362  do_massweight = .false.
1363  if (present(usemasswghtinterp)) do_massweight = usemasswghtinterp
1364 
1365  sv_scale = us%R_to_kg_m3
1366  rl2_t2_to_pa = us%RL2_T2_to_Pa
1367  alpha_ref_mks = alpha_ref * us%kg_m3_to_R
1368 
1369  do n = 1, 5 ! Note that these are reversed from int_density_dz.
1370  wt_t(n) = 0.25 * real(n-1)
1371  wt_b(n) = 1.0 - wt_t(n)
1372  enddo
1373 
1374  ! 1. Compute vertical integrals
1375  do j=jsq,jeq+1; do i=isq,ieq+1
1376  dp = p_b(i,j) - p_t(i,j)
1377  do n=1,5 ! T, S and p are linearly interpolated in the vertical.
1378  p5(n) = rl2_t2_to_pa * (wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j))
1379  s5(n) = wt_t(n) * s_t(i,j) + wt_b(n) * s_b(i,j)
1380  t5(n) = wt_t(n) * t_t(i,j) + wt_b(n) * t_b(i,j)
1381  enddo
1382  if (sv_scale /= 1.0) then
1383  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks, scale=sv_scale)
1384  else
1385  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref_mks)
1386  endif
1387 
1388  ! Use Boole's rule to estimate the interface height anomaly change.
1389  alpha_anom = c1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3))
1390  dza(i,j) = dp*alpha_anom
1391  ! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
1392  ! the interface height anomaly.
1393  if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
1394  (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
1395  enddo ; enddo
1396 
1397  ! 2. Compute horizontal integrals in the x direction
1398  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
1399  ! hWght is the distance measure by which the cell is violation of
1400  ! hydrostatic consistency. For large hWght we bias the interpolation
1401  ! of T,S along the top and bottom integrals, almost like thickness
1402  ! weighting. Note: To work in terrain following coordinates we could
1403  ! offset this distance by the layer thickness to replicate other models.
1404  hwght = 0.0
1405  if (do_massweight) &
1406  hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
1407  if (hwght > 0.) then
1408  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
1409  hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
1410  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1411  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1412  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1413  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1414  else
1415  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1416  endif
1417 
1418  do m=2,4
1419  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1420  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1421 
1422  ! T, S, and p are interpolated in the horizontal. The p interpolation
1423  ! is linear, but for T and S it may be thickness weighted.
1424  p_top = wt_l*p_t(i,j) + wt_r*p_t(i+1,j)
1425  p_bot = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
1426  t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i+1,j)
1427  t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i+1,j)
1428  s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i+1,j)
1429  s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i+1,j)
1430  dp_90(m) = c1_90*(p_bot - p_top)
1431 
1432  ! Salinity, temperature and pressure with linear interpolation in the vertical.
1433  pos = (m-2)*5
1434  do n=1,5
1435  p15(pos+n) = rl2_t2_to_pa * (wt_t(n) * p_top + wt_b(n) * p_bot)
1436  s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
1437  t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
1438  enddo
1439  enddo
1440 
1441  if (sv_scale /= 1.0) then
1442  call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks, scale=sv_scale)
1443  else
1444  call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks)
1445  endif
1446 
1447  intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
1448  do m=2,4
1449  ! Use Boole's rule to estimate the interface height anomaly change.
1450  ! The integrals at the ends of the segment are already known.
1451  pos = (m-2)*5
1452  intp(m) = dp_90(m)*((7.0*(a15(pos+1)+a15(pos+5)) + &
1453  32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
1454  enddo
1455  ! Use Boole's rule to integrate the interface height anomaly values in x.
1456  intx_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
1457  12.0*intp(3))
1458  enddo ; enddo ; endif
1459 
1460  ! 3. Compute horizontal integrals in the y direction
1461  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
1462  ! hWght is the distance measure by which the cell is violation of
1463  ! hydrostatic consistency. For large hWght we bias the interpolation
1464  ! of T,S along the top and bottom integrals, like thickness weighting.
1465  hwght = 0.0
1466  if (do_massweight) &
1467  hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
1468  if (hwght > 0.) then
1469  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
1470  hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
1471  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1472  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1473  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1474  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1475  else
1476  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1477  endif
1478 
1479  do m=2,4
1480  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1481  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1482 
1483  ! T, S, and p are interpolated in the horizontal. The p interpolation
1484  ! is linear, but for T and S it may be thickness weighted.
1485  p_top = wt_l*p_t(i,j) + wt_r*p_t(i,j+1)
1486  p_bot = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
1487  t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i,j+1)
1488  t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i,j+1)
1489  s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i,j+1)
1490  s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i,j+1)
1491  dp_90(m) = c1_90*(p_bot - p_top)
1492 
1493  ! Salinity, temperature and pressure with linear interpolation in the vertical.
1494  pos = (m-2)*5
1495  do n=1,5
1496  p15(pos+n) = rl2_t2_to_pa * (wt_t(n) * p_top + wt_b(n) * p_bot)
1497  s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
1498  t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
1499  enddo
1500  enddo
1501 
1502  if (sv_scale /= 1.0) then
1503  call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks, scale=sv_scale)
1504  else
1505  call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref_mks)
1506  endif
1507 
1508  intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
1509  do m=2,4
1510  ! Use Boole's rule to estimate the interface height anomaly change.
1511  ! The integrals at the ends of the segment are already known.
1512  pos = (m-2)*5
1513  intp(m) = dp_90(m) * ((7.0*(a15(pos+1)+a15(pos+5)) + &
1514  32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
1515  enddo
1516  ! Use Boole's rule to integrate the interface height anomaly values in x.
1517  inty_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
1518  12.0*intp(3))
1519  enddo ; enddo ; endif
1520 
1521 end subroutine int_spec_vol_dp_generic_plm
1522 
1523 
1524 !> Find the depth at which the reconstructed pressure matches P_tgt
1525 subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, &
1526  rho_ref, G_e, EOS, US, P_b, z_out, z_tol)
1527  real, intent(in) :: t_t !< Potential temperature at the cell top [degC]
1528  real, intent(in) :: t_b !< Potential temperature at the cell bottom [degC]
1529  real, intent(in) :: s_t !< Salinity at the cell top [ppt]
1530  real, intent(in) :: s_b !< Salinity at the cell bottom [ppt]
1531  real, intent(in) :: z_t !< Absolute height of top of cell [Z ~> m] (Boussinesq ????)
1532  real, intent(in) :: z_b !< Absolute height of bottom of cell [Z ~> m]
1533  real, intent(in) :: p_t !< Anomalous pressure of top of cell, relative
1534  !! to g*rho_ref*z_t [R L2 T-2 ~> Pa]
1535  real, intent(in) :: p_tgt !< Target pressure at height z_out, relative
1536  !! to g*rho_ref*z_out [R L2 T-2 ~> Pa]
1537  real, intent(in) :: rho_ref !< Reference density with which calculation
1538  !! are anomalous to [R ~> kg m-3]
1539  real, intent(in) :: g_e !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
1540  type(eos_type), pointer :: eos !< Equation of state structure
1541  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1542  real, intent(out) :: p_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa]
1543  real, intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m]
1544  real, optional, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m]
1545 
1546  ! Local variables
1547  real :: dp ! Pressure thickness of the layer [R L2 T-2 ~> Pa]
1548  real :: f_guess, f_l, f_r ! Fractional positions [nondim]
1549  real :: gxrho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]
1550  real :: pa, pa_left, pa_right, pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz [R L2 T-2 ~> Pa]
1551  character(len=240) :: msg
1552 
1553  gxrho = g_e * rho_ref
1554 
1555  ! Anomalous pressure difference across whole cell
1556  dp = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, 1.0, eos)
1557 
1558  p_b = p_t + dp ! Anomalous pressure at bottom of cell
1559 
1560  if (p_tgt <= p_t ) then
1561  z_out = z_t
1562  return
1563  endif
1564 
1565  if (p_tgt >= p_b) then
1566  z_out = z_b
1567  return
1568  endif
1569 
1570  f_l = 0.
1571  pa_left = p_t - p_tgt ! Pa_left < 0
1572  f_r = 1.
1573  pa_right = p_b - p_tgt ! Pa_right > 0
1574  pa_tol = gxrho * 1.0e-5*us%m_to_Z
1575  if (present(z_tol)) pa_tol = gxrho * z_tol
1576 
1577  f_guess = f_l - pa_left / (pa_right - pa_left) * (f_r - f_l)
1578  pa = pa_right - pa_left ! To get into iterative loop
1579  do while ( abs(pa) > pa_tol )
1580 
1581  z_out = z_t + ( z_b - z_t ) * f_guess
1582  pa = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, f_guess, eos) - ( p_tgt - p_t )
1583 
1584  if (pa<pa_left) then
1585  write(msg,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
1586  call mom_error(fatal, 'find_depth_of_pressure_in_cell out of bounds negative: /n'//msg)
1587  elseif (pa<0.) then
1588  pa_left = pa
1589  f_l = f_guess
1590  elseif (pa>pa_right) then
1591  write(msg,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
1592  call mom_error(fatal, 'find_depth_of_pressure_in_cell out of bounds positive: /n'//msg)
1593  elseif (pa>0.) then
1594  pa_right = pa
1595  f_r = f_guess
1596  else ! Pa == 0
1597  return
1598  endif
1599  f_guess = f_l - pa_left / (pa_right - pa_left) * (f_r - f_l)
1600 
1601  enddo
1602 
1603 end subroutine find_depth_of_pressure_in_cell
1604 
1605 
1606 !> Returns change in anomalous pressure change from top to non-dimensional
1607 !! position pos between z_t and z_b
1608 real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS)
1609  real, intent(in) :: t_t !< Potential temperature at the cell top [degC]
1610  real, intent(in) :: t_b !< Potential temperature at the cell bottom [degC]
1611  real, intent(in) :: s_t !< Salinity at the cell top [ppt]
1612  real, intent(in) :: s_b !< Salinity at the cell bottom [ppt]
1613  real, intent(in) :: z_t !< The geometric height at the top of the layer [Z ~> m]
1614  real, intent(in) :: z_b !< The geometric height at the bottom of the layer [Z ~> m]
1615  real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted out to
1616  !! reduce the magnitude of each of the integrals.
1617  real, intent(in) :: g_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
1618  real, intent(in) :: pos !< The fractional vertical position, 0 to 1 [nondim]
1619  type(eos_type), pointer :: eos !< Equation of state structure
1620  real :: fract_dp_at_pos !< The change in pressure from the layer top to
1621  !! fractional position pos [R L2 T-2 ~> Pa]
1622  ! Local variables
1623  real, parameter :: c1_90 = 1.0/90.0 ! A rational constant [nondim]
1624  real :: dz ! Distance from the layer top [Z ~> m]
1625  real :: top_weight, bottom_weight ! Fractional weights at quadrature points [nondim]
1626  real :: rho_ave ! Average density [R ~> kg m-3]
1627  real, dimension(5) :: t5 ! Temperatures at quadrature points [degC]
1628  real, dimension(5) :: s5 ! Salinities at quadrature points [ppt]
1629  real, dimension(5) :: p5 ! Pressures at quadrature points [R L2 T-2 ~> Pa]
1630  real, dimension(5) :: rho5 ! Densities at quadrature points [R ~> kg m-3]
1631  integer :: n
1632 
1633  do n=1,5
1634  ! Evaluate density at five quadrature points
1635  bottom_weight = 0.25*real(n-1) * pos
1636  top_weight = 1.0 - bottom_weight
1637  ! Salinity and temperature points are linearly interpolated
1638  s5(n) = top_weight * s_t + bottom_weight * s_b
1639  t5(n) = top_weight * t_t + bottom_weight * t_b
1640  p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( g_e * rho_ref )
1641  enddo
1642  call calculate_density(t5, s5, p5, rho5, eos)
1643  rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref
1644 
1645  ! Use Boole's rule to estimate the average density
1646  rho_ave = c1_90*(7.0*(rho5(1)+rho5(5)) + 32.0*(rho5(2)+rho5(4)) + 12.0*rho5(3))
1647 
1648  dz = ( z_t - z_b ) * pos
1649  frac_dp_at_pos = g_e * dz * rho_ave
1650 end function frac_dp_at_pos
1651 
1652 end module mom_density_integrals
1653 
1654 !> \namespace mom_density_integrals
1655 !!
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:80
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
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_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:108
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:16
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_density_integrals
Provides integrals of density.
Definition: MOM_density_integrals.F90:2
mom_eos::calculate_specific_vol_derivs
Calculate the derivatives of specific volume with temperature and salinity from T,...
Definition: MOM_EOS.F90:87
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_eos::calculate_spec_vol
Calculates specific volume of sea water from T, S and P.
Definition: MOM_EOS.F90:75
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:68