src/core/MOM_density_integrals.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> Provides integrals of density
6		module MOM_density_integrals
7
8		use MOM_EOS, only : EOS_type
9		use MOM_EOS, only : EOS_quadrature, EOS_domain
10		use MOM_EOS, only : analytic_int_density_dz
11		use MOM_EOS, only : analytic_int_specific_vol_dp
12		use MOM_EOS, only : calculate_density
13		use MOM_EOS, only : calculate_spec_vol
14		use MOM_EOS, only : calculate_specific_vol_derivs
15		use MOM_EOS, only : average_specific_vol
16		use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg
17		use MOM_hor_index, only : hor_index_type
18		use MOM_string_functions, only : uppercase
19		use MOM_variables, only : thermo_var_ptrs
20		use MOM_unit_scaling, only : unit_scale_type
21		use MOM_verticalGrid, only : verticalGrid_type
22
23		implicit none ; private
24
25		#include <MOM_memory.h>
26
27		public int_density_dz
28		public int_density_dz_generic_pcm
29		public int_density_dz_generic_plm
30		public int_density_dz_generic_ppm
31		public int_specific_vol_dp
32		public int_spec_vol_dp_generic_pcm
33		public int_spec_vol_dp_generic_plm
34		public avg_specific_vol
35		public find_depth_of_pressure_in_cell
36		public diagnose_mass_weight_Z, diagnose_mass_weight_p
37
38		contains
39
40		!> Calls the appropriate subroutine to calculate analytical and nearly-analytical
41		!! integrals in z across layers of pressure anomalies, which are
42		!! required for calculating the finite-volume form pressure accelerations in a
43		!! Boussinesq model.
44	0	subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, &
45	0	intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p, &
46		MassWghtInterpVanOnly, h_nv)
47		type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the arrays
48		real, dimension(SZI_(HI),SZJ_(HI)), &
49		intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
50		real, dimension(SZI_(HI),SZJ_(HI)), &
51		intent(in) :: S !< Salinity [S ~> ppt]
52		real, dimension(SZI_(HI),SZJ_(HI)), &
53		intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]
54		real, dimension(SZI_(HI),SZJ_(HI)), &
55		intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]
56		real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is
57		!! subtracted out to reduce the magnitude of each of the
58		!! integrals.
59		real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used
60		!! to calculate the pressure (as p~=-zrho_0G_e)
61		!! used in the equation of state.
62		real, intent(in) :: G_e !< The Earth's gravitational acceleration
63		!! [L2 Z-1 T-2 ~> m s-2]
64		type(EOS_type), intent(in) :: EOS !< Equation of state structure
65		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
66		real, dimension(SZI_(HI),SZJ_(HI)), &
67		intent(inout) :: dpa !< The change in the pressure anomaly
68		!! across the layer [R L2 T-2 ~> Pa]
69		real, dimension(SZI_(HI),SZJ_(HI)), &
70		optional, intent(inout) :: intz_dpa !< The integral through the thickness of the
71		!! layer of the pressure anomaly relative to the
72		!! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]
73		real, dimension(SZIB_(HI),SZJ_(HI)), &
74		optional, intent(inout) :: intx_dpa !< The integral in x of the difference between
75		!! the pressure anomaly at the top and bottom of the
76		!! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
77		real, dimension(SZI_(HI),SZJB_(HI)), &
78		optional, intent(inout) :: inty_dpa !< The integral in y of the difference between
79		!! the pressure anomaly at the top and bottom of the
80		!! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
81		real, dimension(SZI_(HI),SZJ_(HI)), &
82		optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]
83		real, dimension(SZI_(HI),SZJ_(HI)), &
84		optional, intent(in) :: SSH !< The sea surface height [Z ~> m]
85		real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]
86		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
87		!! mass weighting to interpolate T/S in integrals
88		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
89		!! of T/S unless one side smaller than h_nv (i.e. vanished)
90		real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]
91
92		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
93		optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
94
95	0	if (EOS_quadrature(EOS)) then
96		call int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, &
97		intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, &
98		MassWghtInterp, Z_0p=Z_0p, MassWghtInterpVanOnly=MassWghtInterpVanOnly, &
99	0	h_nv=h_nv)
100		else
101		call analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, &
102	0	intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p=Z_0p)
103		endif
104
105	0	end subroutine int_density_dz
106
107
108		!> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which
109		!! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model.
110	0	subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
111	0	EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &
112	0	dz_neglect, MassWghtInterp, use_inaccurate_form, Z_0p, &
113		MassWghtInterpVanOnly, h_nv)
114		type(hor_index_type), intent(in) :: HI !< Horizontal index type for input variables.
115		real, dimension(SZI_(HI),SZJ_(HI)), &
116		intent(in) :: T !< Potential temperature of the layer [C ~> degC]
117		real, dimension(SZI_(HI),SZJ_(HI)), &
118		intent(in) :: S !< Salinity of the layer [S ~> ppt]
119		real, dimension(SZI_(HI),SZJ_(HI)), &
120		intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]
121		real, dimension(SZI_(HI),SZJ_(HI)), &
122		intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]
123		real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is
124		!! subtracted out to reduce the magnitude
125		!! of each of the integrals.
126		real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used
127		!! to calculate the pressure (as p~=-zrho_0G_e)
128		!! used in the equation of state.
129		real, intent(in) :: G_e !< The Earth's gravitational acceleration
130		!! [L2 Z-1 T-2 ~> m s-2]
131		type(EOS_type), intent(in) :: EOS !< Equation of state structure
132		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
133		real, dimension(SZI_(HI),SZJ_(HI)), &
134		intent(inout) :: dpa !< The change in the pressure anomaly
135		!! across the layer [R L2 T-2 ~> Pa]
136		real, dimension(SZI_(HI),SZJ_(HI)), &
137		optional, intent(inout) :: intz_dpa !< The integral through the thickness of the
138		!! layer of the pressure anomaly relative to the
139		!! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]
140		real, dimension(SZIB_(HI),SZJ_(HI)), &
141		optional, intent(inout) :: intx_dpa !< The integral in x of the difference between
142		!! the pressure anomaly at the top and bottom of the
143		!! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
144		real, dimension(SZI_(HI),SZJB_(HI)), &
145		optional, intent(inout) :: inty_dpa !< The integral in y of the difference between
146		!! the pressure anomaly at the top and bottom of the
147		!! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
148		real, dimension(SZI_(HI),SZJ_(HI)), &
149		optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]
150		real, dimension(SZI_(HI),SZJ_(HI)), &
151		optional, intent(in) :: SSH !< The sea surface height [Z ~> m]
152		real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]
153		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
154		!! mass weighting to interpolate T/S in integrals
155		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
156		!! of T/S unless one side smaller than h_nv (i.e. vanished)
157		real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]
158		logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of
159		!! density anomalies, as was used prior to March 2018.
160		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
161		optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
162
163		! Local variables
164	0	real :: T5((5HI%iscB+1):(5(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC]
165	0	real :: S5((5HI%iscB+1):(5(HI%iecB+2))) ! Salinities along a line of subgrid locations [S ~> ppt]
166	0	real :: p5((5HI%iscB+1):(5(HI%iecB+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]
167	0	real :: r5((5HI%iscB+1):(5(HI%iecB+2))) ! Densities anomalies along a line of subgrid locations [R ~> kg m-3]
168	0	real :: T15((15HI%iscB+1):(15(HI%iecB+1))) ! Temperatures at an array of subgrid locations [C ~> degC]
169	0	real :: S15((15HI%iscB+1):(15(HI%iecB+1))) ! Salinities at an array of subgrid locations [S ~> ppt]
170	0	real :: p15((15HI%iscB+1):(15(HI%iecB+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]
171	0	real :: r15((15HI%iscB+1):(15(HI%iecB+1))) ! Densities at an array of subgrid locations [R ~> kg m-3]
172		real :: rho_anom ! The depth averaged density anomaly [R ~> kg m-3]
173		real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
174		real :: GxRho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]
175		real :: dz ! The layer thickness [Z ~> m]
176	0	real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]
177	0	real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]
178	0	real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
179		real :: hWght ! A pressure-thickness below topography [Z ~> m]
180		real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m]
181		real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2]
182		real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim]
183		real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]
184		real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]
185		real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]
186		real :: intz(5) ! The gravitational acceleration times the integrals of density
187		! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]
188		logical :: do_massWeight ! Indicates whether to do mass weighting near bathymetry
189		logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
190		real :: massWeightNVonlyToggle ! A non-dimensional toggle factor for only using mass weighting
191		! if at least one side vanished (0 or 1) [nondim]
192		real :: h_nonvanished ! nonvanished height [Z ~> m]
193		logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation
194		! of density anomalies.
195		integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state
196		integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state
197		integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state
198		integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m, n, pos
199
200		! These array bounds work for the indexing convention of the input arrays, but
201		! on the computational domain defined for the output arrays.
202	0	Isq = HI%IscB ; Ieq = HI%IecB
203	0	Jsq = HI%JscB ; Jeq = HI%JecB
204	0	is = HI%isc ; ie = HI%iec
205	0	js = HI%jsc ; je = HI%jec
206
207	0	GxRho = G_e * rho_0
208	0	if (present(Z_0p)) then
209	0	do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
210	0	z0pres(i,j) = Z_0p(i,j)
211		enddo ; enddo
212		else
213	0	z0pres(:,:) = 0.0
214		endif
215	0	use_rho_ref = .true.
216	0	if (present(use_inaccurate_form)) then
217	0	if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form
218		endif
219
220	0	do_massWeight = .false. ; top_massWeight = .false.
221	0	if (present(MassWghtInterp)) then
222	0	do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
223	0	top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
224	0	if (do_massWeight .and. .not.present(bathyT)) call MOM_error(FATAL, &
225	0	"int_density_dz_generic: bathyT must be present if near-bottom mass weighting is in use.")
226	0	if (top_massWeight .and. .not.present(SSH)) call MOM_error(FATAL, &
227	0	"int_density_dz_generic: SSH must be present if near-surface mass weighting is in use.")
228	0	if ((do_massWeight .or. top_massWeight) .and. .not.present(dz_neglect)) call MOM_error(FATAL, &
229	0	"int_density_dz_generic: dz_neglect must be present if mass weighting is in use.")
230		endif
231	0	massWeightNVonlyToggle = 1.
232	0	if (present(MassWghtInterpVanOnly)) then
233	0	if (MassWghtInterpVanOnly) massWeightNVonlyToggle = 0.
234		endif
235	0	h_nonvanished = 0.
236	0	if (present(h_nv)) then
237	0	h_nonvanished = h_nv
238		endif
239
240		! Set the loop ranges for equation of state calculations at various points.
241	0	EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2)
242	0	EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1)
243	0	EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1)
244
245	0	do j=Jsq,Jeq+1
246	0	do i=Isq,Ieq+1
247	0	dz = z_t(i,j) - z_b(i,j)
248	0	do n=1,5
249	0	T5(i5+n) = T(i,j) ; S5(i5+n) = S(i,j)
250	0	p5(i5+n) = -GxRho((z_t(i,j) - z0pres(i,j)) - 0.25real(n-1)dz)
251		enddo
252		enddo
253
254	0	if (use_rho_ref) then
255	0	call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5, rho_ref=rho_ref)
256		else
257	0	call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5)
258		endif
259
260	0	do i=Isq,Ieq+1
261		! Use Boole's rule to estimate the pressure anomaly change.
262	0	rho_anom = C1_90(7.0(r5(i5+1)+r5(i5+5)) + 32.0(r5(i5+2)+r5(i5+4)) + 12.0r5(i*5+3))
263	0	if (.not.use_rho_ref) rho_anom = rho_anom - rho_ref
264	0	dz = z_t(i,j) - z_b(i,j)
265	0	dpa(i,j) = G_edzrho_anom
266		! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
267		! the pressure anomaly.
268	0	if (present(intz_dpa)) intz_dpa(i,j) = 0.5G_edz*2 &
269	0	(rho_anom - C1_90(16.0(r5(i5+4)-r5(i5+2)) + 7.0(r5(i5+5)-r5(i*5+1))) )
270		enddo
271		enddo
272
273	0	if (present(intx_dpa)) then ; do j=js,je
274	0	do I=Isq,Ieq
275		! hWght is the distance measure by which the cell is violation of
276		! hydrostatic consistency. For large hWght we bias the interpolation of
277		! T & S along the top and bottom integrals, akin to thickness weighting.
278	0	hWght = 0.0
279	0	if (do_massWeight) &
280	0	hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j))
281	0	if (top_massWeight) &
282	0	hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j))
283		! If both sides are nonvanished, then set it back to zero.
284	0	if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i+1,j) - z_b(i+1,j)) > h_nonvanished)) then
285	0	hWght = massWeightNVonlyToggle * hWght
286		endif
287	0	if (hWght > 0.) then
288	0	hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
289	0	hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
290	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
291	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
292	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
293	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
294		else
295	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
296		endif
297
298	0	do m=2,4
299		! T, S, and z are interpolated in the horizontal. The z interpolation
300		! is linear, but for T and S it may be thickness weighted.
301	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
302	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
303	0	dz_x(m,i) = (wt_L(z_t(i,j) - z_b(i,j))) + (wt_R(z_t(i+1,j) - z_b(i+1,j)))
304	0	pos = i15+(m-2)5
305	0	T15(pos+1) = (wtT_LT(i,j)) + (wtT_RT(i+1,j))
306	0	S15(pos+1) = (wtT_LS(i,j)) + (wtT_RS(i+1,j))
307	0	p15(pos+1) = -GxRho * ((wt_L(z_t(i,j)-z0pres(i,j))) + (wt_R(z_t(i+1,j)-z0pres(i+1,j))))
308	0	do n=2,5
309	0	T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
310	0	p15(pos+n) = p15(pos+n-1) + GxRho0.25dz_x(m,i)
311		enddo
312		enddo
313		enddo
314
315	0	if (use_rho_ref) then
316	0	call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15, rho_ref=rho_ref)
317		else
318	0	call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15)
319		endif
320
321	0	do I=Isq,Ieq
322	0	intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
323		! Use Boole's rule to estimate the pressure anomaly change.
324	0	if (use_rho_ref) then
325	0	do m=2,4
326	0	pos = i15+(m-2)5
327		intz(m) = (G_edz_x(m,i)(C1_90( 7.0(r15(pos+1)+r15(pos+5)) + &
328		32.0*(r15(pos+2)+r15(pos+4)) + &
329	0	12.0*r15(pos+3)) ))
330		enddo
331		else
332	0	do m=2,4
333	0	pos = i15+(m-2)5
334		intz(m) = (G_edz_x(m,i)(C1_90( 7.0(r15(pos+1)+r15(pos+5)) + &
335		32.0*(r15(pos+2)+r15(pos+4)) + &
336	0	12.0*r15(pos+3)) - rho_ref ))
337		enddo
338		endif
339		! Use Boole's rule to integrate the bottom pressure anomaly values in x.
340		intx_dpa(i,j) = C1_90(7.0(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
341	0	12.0*intz(3))
342		enddo
343		enddo ; endif
344
345	0	if (present(inty_dpa)) then ; do J=Jsq,Jeq
346	0	do i=is,ie
347		! hWght is the distance measure by which the cell is violation of
348		! hydrostatic consistency. For large hWght we bias the interpolation of
349		! T & S along the top and bottom integrals, akin to thickness weighting.
350	0	hWght = 0.0
351	0	if (do_massWeight) &
352	0	hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j))
353	0	if (top_massWeight) &
354	0	hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1))
355		! If both sides are nonvanished, then set it back to zero.
356	0	if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i,j+1) - z_b(i,j+1)) > h_nonvanished)) then
357	0	hWght = massWeightNVonlyToggle * hWght
358		endif
359	0	if (hWght > 0.) then
360	0	hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
361	0	hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
362	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
363	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
364	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
365	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
366		else
367	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
368		endif
369
370	0	do m=2,4
371		! T, S, and z are interpolated in the horizontal. The z interpolation
372		! is linear, but for T and S it may be thickness weighted.
373	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
374	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
375	0	dz_y(m,i) = (wt_L(z_t(i,j) - z_b(i,j))) + (wt_R(z_t(i,j+1) - z_b(i,j+1)))
376	0	pos = i15+(m-2)5
377	0	T15(pos+1) = (wtT_LT(i,j)) + (wtT_RT(i,j+1))
378	0	S15(pos+1) = (wtT_LS(i,j)) + (wtT_RS(i,j+1))
379	0	p15(pos+1) = -GxRho * ((wt_L(z_t(i,j)-z0pres(i,j))) + (wt_R(z_t(i,j+1)-z0pres(i,j+1))))
380	0	do n=2,5
381	0	T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
382	0	p15(pos+n) = p15(pos+n-1) + GxRho0.25dz_y(m,i)
383		enddo
384		enddo
385		enddo
386
387	0	if (use_rho_ref) then
388		call calculate_density(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
389	0	r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref)
390		else
391		call calculate_density(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
392	0	r15(15*HI%isc+1:), EOS, EOSdom_h15)
393		endif
394
395	0	do i=is,ie
396	0	intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
397		! Use Boole's rule to estimate the pressure anomaly change.
398	0	do m=2,4
399	0	pos = i15+(m-2)5
400	0	if (use_rho_ref) then
401		intz(m) = (G_edz_y(m,i)(C1_90(7.0(r15(pos+1)+r15(pos+5)) + &
402		32.0*(r15(pos+2)+r15(pos+4)) + &
403	0	12.0*r15(pos+3)) ))
404		else
405		intz(m) = (G_edz_y(m,i)(C1_90(7.0(r15(pos+1)+r15(pos+5)) + &
406		32.0*(r15(pos+2)+r15(pos+4)) + &
407	0	12.0*r15(pos+3)) - rho_ref ))
408		endif
409		enddo
410		! Use Boole's rule to integrate the values.
411		inty_dpa(i,j) = C1_90(7.0(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
412	0	12.0*intz(3))
413		enddo
414		enddo ; endif
415	0	end subroutine int_density_dz_generic_pcm
416
417
418		!> Compute pressure gradient force integrals by quadrature for the case where
419		!! T and S are linear profiles.
420	3600	subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
421	3600	rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, dpa, &
422	5400	intz_dpa, intx_dpa, inty_dpa, MassWghtInterp, &
423	1800	use_inaccurate_form, Z_0p, MassWghtInterpVanOnly, h_nv)
424		integer, intent(in) :: k !< Layer index to calculate integrals for
425		type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays
426		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
427		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
428		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
429		intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC]
430		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
431		intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC]
432		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
433		intent(in) :: S_t !< Salinity at the cell top [S ~> ppt]
434		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
435		intent(in) :: S_b !< Salinity at the cell bottom [S ~> ppt]
436		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &
437		intent(in) :: e !< Height of interfaces [Z ~> m]
438		real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted
439		!! out to reduce the magnitude of each of the integrals.
440		real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used to calculate
441		!! the pressure (as p~=-zrho_0G_e) used in the equation of state.
442		real, intent(in) :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
443		real, intent(in) :: dz_subroundoff !< A minuscule thickness change [Z ~> m]
444		real, dimension(SZI_(HI),SZJ_(HI)), &
445		intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]
446		type(EOS_type), intent(in) :: EOS !< Equation of state structure
447		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
448		logical, intent(in) :: use_stanley_eos !< If true, turn on Stanley SGS T variance parameterization
449		real, dimension(SZI_(HI),SZJ_(HI)), &
450		intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]
451		real, dimension(SZI_(HI),SZJ_(HI)), &
452		optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of
453		!! the pressure anomaly relative to the anomaly at the
454		!! top of the layer [R L2 Z T-2 ~> Pa m]
455		real, dimension(SZIB_(HI),SZJ_(HI)), &
456		optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
457		!! pressure anomaly at the top and bottom of the layer
458		!! divided by the x grid spacing [R L2 T-2 ~> Pa]
459		real, dimension(SZI_(HI),SZJB_(HI)), &
460		optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
461		!! pressure anomaly at the top and bottom of the layer
462		!! divided by the y grid spacing [R L2 T-2 ~> Pa]
463		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
464		!! mass weighting to interpolate T/S in integrals
465		logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of
466		!! density anomalies, as was used prior to March 2018.
467		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
468		!! of T/S unless one side smaller than h_nv (i.e. vanished)
469		real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]
470		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
471		optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
472
473		! This subroutine calculates (by numerical quadrature) integrals of
474		! pressure anomalies across layers, which are required for calculating the
475		! finite-volume form pressure accelerations in a Boussinesq model. The one
476		! potentially dodgy assumption here is that rho_0 is used both in the denominator
477		! of the accelerations, and in the pressure used to calculated density (the
478		! latter being -zrho_0G_e). These two uses could be separated if need be.
479		!
480		! It is assumed that the salinity and temperature profiles are linear in the
481		! vertical. The top and bottom values within each layer are provided and
482		! a linear interpolation is used to compute intermediate values.
483
484		! Local variables
485	3600	real :: T5((5HI%iscB+1):(5(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC]
486	3600	real :: S5((5HI%iscB+1):(5(HI%iecB+2))) ! Salinities along a line of subgrid locations [S ~> ppt]
487	3600	real :: T25((5HI%iscB+1):(5(HI%iecB+2))) ! SGS temperature variance along a line of subgrid
488		! locations [C2 ~> degC2]
489	3600	real :: TS5((5HI%iscB+1):(5(HI%iecB+2))) ! SGS temp-salt covariance along a line of subgrid
490		! locations [C S ~> degC ppt]
491	3600	real :: S25((5HI%iscB+1):(5(HI%iecB+2))) ! SGS salinity variance along a line of subgrid locations [S2 ~> ppt2]
492	3600	real :: p5((5HI%iscB+1):(5(HI%iecB+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]
493	3600	real :: r5((5HI%iscB+1):(5(HI%iecB+2))) ! Densities anomalies along a line of subgrid
494		! locations [R ~> kg m-3]
495	3600	real :: u5((5HI%iscB+1):(5(HI%iecB+2))) ! Densities anomalies along a line of subgrid locations
496		! (used for inaccurate form) [R ~> kg m-3]
497	3600	real :: T15((15HI%iscB+1):(15(HI%iecB+1))) ! Temperatures at an array of subgrid locations [C ~> degC]
498	3600	real :: S15((15HI%iscB+1):(15(HI%iecB+1))) ! Salinities at an array of subgrid locations [S ~> ppt]
499	3600	real :: T215((15HI%iscB+1):(15(HI%iecB+1))) ! SGS temperature variance along a line of subgrid
500		! locations [C2 ~> degC2]
501	3600	real :: TS15((15HI%iscB+1):(15(HI%iecB+1))) ! SGS temp-salt covariance along a line of subgrid
502		! locations [C S ~> degC ppt]
503	3600	real :: S215((15HI%iscB+1):(15(HI%iecB+1))) ! SGS salinity variance along a line of subgrid
504		! locations [S2 ~> ppt2]
505	3600	real :: p15((15HI%iscB+1):(15(HI%iecB+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]
506	3600	real :: r15((15HI%iscB+1):(15(HI%iecB+1))) ! Densities at an array of subgrid locations [R ~> kg m-3]
507		real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]
508		real :: rho_anom ! A density anomaly [R ~> kg m-3]
509		real :: w_left, w_right ! Left and right weights [nondim]
510		real :: intz(5) ! The gravitational acceleration times the integrals of density
511		! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]
512		real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
513		real :: GxRho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]
514	3600	real :: dz(HI%iscB:HI%iecB+1) ! Layer thicknesses at tracer points [Z ~> m]
515	3600	real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]
516	3600	real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]
517		real :: massWeightToggle ! A non-dimensional toggle factor for near-bottom mass weighting (0 or 1) [nondim]
518		real :: TopWeightToggle ! A non-dimensional toggle factor for near-surface mass weighting (0 or 1) [nondim]
519		real :: massWeightNVonlyToggle ! A non-dimensional toggle factor for only using mass weighting
520		! if at least one side vanished (0 or 1) [nondim]
521		real :: Ttl, Tbl, Ttr, Tbr ! Temperatures at the velocity cell corners [C ~> degC]
522		real :: Stl, Sbl, Str, Sbr ! Salinities at the velocity cell corners [S ~> ppt]
523	3600	real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
524		real :: hWght ! A topographically limited thickness weight [Z ~> m]
525		real :: hWghtTop ! An ice draft limited thickness weight [Z ~> m]
526		real :: hL, hR ! Thicknesses to the left and right [Z ~> m]
527		real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
528		real :: h_nonvanished ! nonvanished height [Z ~> m]
529		logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation
530		! of density anomalies.
531		logical :: use_varT, use_varS, use_covarTS ! Logicals for SGS variances fields
532		integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state
533		integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state
534		integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state
535		integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos
536
537	1800	Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB
538
539	1800	GxRho = G_e * rho_0
540	1800	if (present(Z_0p)) then
541	13728600	do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
542	13726800	z0pres(i,j) = Z_0p(i,j)
543		enddo ; enddo
544		else
545	0	z0pres(:,:) = 0.0
546		endif
547	1800	massWeightToggle = 0. ; TopWeightToggle = 0.
548	1800	if (present(MassWghtInterp)) then
549	1800	if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1.
550	1800	if (BTEST(MassWghtInterp, 1)) TopWeightToggle = 1.
551		endif
552	1800	massWeightNVonlyToggle = 1.
553	1800	if (present(MassWghtInterpVanOnly)) then
554	1800	if (MassWghtInterpVanOnly) massWeightNVonlyToggle = 0.
555		endif
556	1800	h_nonvanished = 0.
557	1800	if (present(h_nv)) then
558	1800	h_nonvanished = h_nv
559		endif
560	1800	use_rho_ref = .true.
561	1800	if (present(use_inaccurate_form)) use_rho_ref = .not. use_inaccurate_form
562
563	1800	use_varT = .false. !ensure initialized
564	1800	use_covarTS = .false.
565	1800	use_varS = .false.
566	1800	if (use_stanley_eos) then
567	0	use_varT = associated(tv%varT)
568	0	use_covarTS = associated(tv%covarTS)
569	0	use_varS = associated(tv%varS)
570		endif
571
572	1099800	T25(:) = 0.
573	1099800	TS5(:) = 0.
574	1099800	S25(:) = 0.
575	3268800	T215(:) = 0.
576	3268800	TS15(:) = 0.
577	3268800	S215(:) = 0.
578
579	10800	do n = 1, 5
580	9000	wt_t(n) = 0.25 * real(5-n)
581	10800	wt_b(n) = 1.0 - wt_t(n)
582		enddo
583
584		! Set the loop ranges for equation of state calculations at various points.
585	1800	EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2)
586	1800	EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1)
587	1800	EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1)
588
589		! 1. Compute vertical integrals
590	113400	do j=Jsq,Jeq+1
591	13726800	do i = Isq,Ieq+1
592	13615200	dz(i) = e(i,j,K) - e(i,j,K+1)
593	81691200	do n=1,5
594	68076000	p5(i5+n) = -GxRho((e(i,j,K) - z0pres(i,j)) - 0.25real(n-1)dz(i))
595		! Salinity and temperature points are linearly interpolated
596	68076000	S5(i5+n) = wt_t(n) S_t(i,j,k) + wt_b(n) * S_b(i,j,k)
597	81691200	T5(i5+n) = wt_t(n) T_t(i,j,k) + wt_b(n) * T_b(i,j,k)
598		enddo
599	13615200	if (use_varT) T25(i5+1:i5+5) = tv%varT(i,j,k)
600	13615200	if (use_covarTS) TS5(i5+1:i5+5) = tv%covarTS(i,j,k)
601	13726800	if (use_varS) S25(i5+1:i5+5) = tv%varS(i,j,k)
602		enddo
603	111600	if (use_Stanley_eos) then
604	0	call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, EOSdom_h5, rho_ref=rho_ref)
605		else
606	111600	if (use_rho_ref) then
607	111600	call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5, rho_ref=rho_ref)
608		else
609	0	call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5)
610	0	u5(:) = r5(:) - rho_ref
611		endif
612		endif
613
614	113400	if (use_rho_ref) then
615	13726800	do i=Isq,Ieq+1
616		! Use Boole's rule to estimate the pressure anomaly change.
617	13615200	rho_anom = C1_90(7.0(r5(i5+1)+r5(i5+5)) + 32.0(r5(i5+2)+r5(i5+4)) + 12.0r5(i*5+3))
618	13615200	dpa(i,j) = G_edz(i)rho_anom
619	13726800	if (present(intz_dpa)) then
620		! Use a Boole's-rule-like fifth-order accurate estimate of
621		! the double integral of the pressure anomaly.
622		intz_dpa(i,j) = 0.5G_edz(i)*2 &
623	13615200	(rho_anom - C1_90(16.0(r5(i5+4)-r5(i5+2)) + 7.0(r5(i5+5)-r5(i*5+1))) )
624		endif
625		enddo
626		else
627	0	do i=Isq,Ieq+1
628		! Use Boole's rule to estimate the pressure anomaly change.
629		rho_anom = C1_90(7.0(r5(i5+1)+r5(i5+5)) + 32.0(r5(i5+2)+r5(i5+4)) + 12.0r5(i*5+3)) &
630	0	- rho_ref
631	0	dpa(i,j) = G_edz(i)rho_anom
632	0	if (present(intz_dpa)) then
633		! Use a Boole's-rule-like fifth-order accurate estimate of
634		! the double integral of the pressure anomaly.
635		intz_dpa(i,j) = 0.5G_edz(i)*2 &
636	0	(rho_anom - C1_90(16.0(u5(i5+4)-u5(i5+2)) + 7.0(u5(i5+5)-u5(i*5+1))) )
637		endif
638		enddo
639		endif
640		enddo ! end loops on j
641
642		! 2. Compute horizontal integrals in the x direction
643	109800	if (present(intx_dpa)) then ; do j=HI%jsc,HI%jec
644	13176000	do I=Isq,Ieq
645		! Corner values of T and S
646		! hWght is the distance measure by which the cell is violation of
647		! hydrostatic consistency. For large hWght we bias the interpolation
648		! of T,S along the top and bottom integrals, almost like thickness
649		! weighting.
650		! Note: To work in terrain following coordinates we could offset
651		! this distance by the layer thickness to replicate other models.
652		hWght = massWeightToggle * &
653	13068000	max(0., -bathyT(i,j)-e(i+1,j,K), -bathyT(i+1,j)-e(i,j,K))
654		! CY: The below code just uses top interface, which may be bad in high res open ocean
655		! We want something like if (pa(i+1,k+1)<pa(i,1)) or (pa(i+1,1) <pa(i,k+1)) then...
656		! but pressures are not passed through to this submodule, and tv just has surface press.
657		!if ((p(i+1,j,k+1)<p(i,j,1)).or.(tv%p(i+1,j,k+1)<tv%p(i,j,1))) then
658		hWghtTop = TopWeightToggle * &
659	13068000	max(0., e(i+1,j,K+1)-e(i,j,1), e(i,j,K+1)-e(i+1,j,1))
660		!else ! pressure criteria not activated
661		! hWghtTop = 0.
662		!endif
663		! Set it to be max of the bottom and top hWghts:
664	13068000	hWght = max(hWght, hWghtTop)
665		! If both sides are nonvanished, then set it back to zero.
666	13068000	if (((e(i,j,K) - e(i,j,K+1)) > h_nonvanished) .and. ((e(i+1,j,K) - e(i+1,j,K+1)) > h_nonvanished)) then
667	13068000	hWght = massWeightNVonlyToggle * hWght
668		endif
669	13068000	if (hWght > 0.) then
670	3069590	hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff
671	3069590	hR = (e(i+1,j,K) - e(i+1,j,K+1)) + dz_subroundoff
672	3069590	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
673	3069590	iDenom = 1./( hWght(hR + hL) + hLhR )
674	3069590	Ttl = ( (hWghthR)T_t(i+1,j,k) + (hWghthL + hRhL)T_t(i,j,k) ) iDenom
675	3069590	Ttr = ( (hWghthL)T_t(i,j,k) + (hWghthR + hRhL)T_t(i+1,j,k) ) iDenom
676	3069590	Tbl = ( (hWghthR)T_b(i+1,j,k) + (hWghthL + hRhL)T_b(i,j,k) ) iDenom
677	3069590	Tbr = ( (hWghthL)T_b(i,j,k) + (hWghthR + hRhL)T_b(i+1,j,k) ) iDenom
678	3069590	Stl = ( (hWghthR)S_t(i+1,j,k) + (hWghthL + hRhL)S_t(i,j,k) ) iDenom
679	3069590	Str = ( (hWghthL)S_t(i,j,k) + (hWghthR + hRhL)S_t(i+1,j,k) ) iDenom
680	3069590	Sbl = ( (hWghthR)S_b(i+1,j,k) + (hWghthL + hRhL)S_b(i,j,k) ) iDenom
681	3069590	Sbr = ( (hWghthL)S_b(i,j,k) + (hWghthR + hRhL)S_b(i+1,j,k) ) iDenom
682		else
683	9998410	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)
684	9998410	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)
685		endif
686
687	52380000	do m=2,4
688	39204000	w_left = wt_t(m) ; w_right = wt_b(m)
689	39204000	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)))
690
691		! Salinity and temperature points are linearly interpolated in
692		! the horizontal. The subscript (1) refers to the top value in
693		! the vertical profile while subscript (5) refers to the bottom
694		! value in the vertical profile.
695	39204000	pos = i15+(m-2)5
696	39204000	T15(pos+1) = (w_leftTtl) + (w_rightTtr)
697	39204000	T15(pos+5) = (w_leftTbl) + (w_rightTbr)
698
699	39204000	S15(pos+1) = (w_leftStl) + (w_rightStr)
700	39204000	S15(pos+5) = (w_leftSbl) + (w_rightSbr)
701
702	39204000	p15(pos+1) = -GxRho * ((w_left(e(i,j,K)-z0pres(i,j))) + (w_right(e(i+1,j,K)-z0pres(i+1,j))))
703
704		! Pressure
705	196020000	do n=2,5
706	196020000	p15(pos+n) = p15(pos+n-1) + GxRho0.25dz_x(m,i)
707		enddo
708
709		! Salinity and temperature (linear interpolation in the vertical)
710	156816000	do n=2,4
711	117612000	S15(pos+n) = wt_t(n) * S15(pos+1) + wt_b(n) * S15(pos+5)
712	156816000	T15(pos+n) = wt_t(n) * T15(pos+1) + wt_b(n) * T15(pos+5)
713		enddo
714	39204000	if (use_varT) T215(pos+1:pos+5) = (w_lefttv%varT(i,j,k)) + (w_righttv%varT(i+1,j,k))
715	39204000	if (use_covarTS) TS15(pos+1:pos+5) = (w_lefttv%covarTS(i,j,k)) + (w_righttv%covarTS(i+1,j,k))
716	52272000	if (use_varS) S215(pos+1:pos+5) = (w_lefttv%varS(i,j,k)) + (w_righttv%varS(i+1,j,k))
717		enddo
718		enddo
719
720	108000	if (use_stanley_eos) then
721	0	call calculate_density(T15, S15, p15, T215, TS15, S215, r15, EOS, EOSdom_q15, rho_ref=rho_ref)
722		else
723	108000	if (use_rho_ref) then
724	108000	call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15, rho_ref=rho_ref)
725		else
726	0	call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15)
727		endif
728		endif
729
730	13177800	do I=Isq,Ieq
731	13068000	intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
732
733		! Use Boole's rule to estimate the pressure anomaly change.
734	13068000	if (use_rho_ref) then
735	52272000	do m = 2,4
736	39204000	pos = i15+(m-2)5
737		intz(m) = (G_edz_x(m,i)( C1_90(7.0(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &
738	52272000	12.0*r15(pos+3)) ))
739		enddo
740		else
741	0	do m = 2,4
742	0	pos = i15+(m-2)5
743		intz(m) = (G_edz_x(m,i)( C1_90(7.0(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &
744	0	12.0*r15(pos+3)) - rho_ref ))
745		enddo
746		endif
747		! Use Boole's rule to integrate the bottom pressure anomaly values in x.
748		intx_dpa(I,j) = C1_90(7.0(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
749	13176000	12.0*intz(3))
750		enddo
751		enddo ; endif
752
753		! 3. Compute horizontal integrals in the y direction
754	111600	if (present(inty_dpa)) then ; do J=Jsq,Jeq
755	13285800	do i=HI%isc,HI%iec
756		! Corner values of T and S
757		! hWght is the distance measure by which the cell is violation of
758		! hydrostatic consistency. For large hWght we bias the interpolation
759		! of T,S along the top and bottom integrals, almost like thickness
760		! weighting.
761		! Note: To work in terrain following coordinates we could offset
762		! this distance by the layer thickness to replicate other models.
763		hWght = massWeightToggle * &
764	13176000	max(0., -bathyT(i,j)-e(i,j+1,K), -bathyT(i,j+1)-e(i,j,K))
765		! CY: The below code just uses top interface, which may be bad in high res open ocean
766		! We want something like if (pa(j+1,k+1)<pa(j,1)) or (pa(j+1,1) <pa(i,j,k+1)) then...
767		! but pressures are not passed through to this submodule, and tv just has surface press.
768		!if ((p(i,j+1,k+1)<p(i,j,1)).or.(tv%p(i,j+1,k+1)<tv%p(i,j,1))) then
769		hWghtTop = TopWeightToggle * &
770	13176000	max(0., e(i,j+1,K+1)-e(i,j,1), e(i,j,K+1)-e(i,j+1,1))
771		!else ! pressure criteria not activated
772		! hWghtTop = 0.
773		!endif
774		! Set it to be max of the bottom and top hWghts:
775	13176000	hWght = max(hWght, hWghtTop)
776		! If both sides are nonvanished, then set it back to zero.
777	13176000	if (((e(i,j,K) - e(i,j,K+1)) > h_nonvanished) .and. ((e(i,j+1,K) - e(i,j+1,K+1)) > h_nonvanished)) then
778	13176000	hWght = massWeightNVonlyToggle * hWght
779		endif
780
781	13176000	if (hWght > 0.) then
782	3174185	hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff
783	3174185	hR = (e(i,j+1,K) - e(i,j+1,K+1)) + dz_subroundoff
784	3174185	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
785	3174185	iDenom = 1./( hWght(hR + hL) + hLhR )
786	3174185	Ttl = ( (hWghthR)T_t(i,j+1,k) + (hWghthL + hRhL)T_t(i,j,k) ) iDenom
787	3174185	Ttr = ( (hWghthL)T_t(i,j,k) + (hWghthR + hRhL)T_t(i,j+1,k) ) iDenom
788	3174185	Tbl = ( (hWghthR)T_b(i,j+1,k) + (hWghthL + hRhL)T_b(i,j,k) ) iDenom
789	3174185	Tbr = ( (hWghthL)T_b(i,j,k) + (hWghthR + hRhL)T_b(i,j+1,k) ) iDenom
790	3174185	Stl = ( (hWghthR)S_t(i,j+1,k) + (hWghthL + hRhL)S_t(i,j,k) ) iDenom
791	3174185	Str = ( (hWghthL)S_t(i,j,k) + (hWghthR + hRhL)S_t(i,j+1,k) ) iDenom
792	3174185	Sbl = ( (hWghthR)S_b(i,j+1,k) + (hWghthL + hRhL)S_b(i,j,k) ) iDenom
793	3174185	Sbr = ( (hWghthL)S_b(i,j,k) + (hWghthR + hRhL)S_b(i,j+1,k) ) iDenom
794		else
795	10001815	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)
796	10001815	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)
797		endif
798
799	52813800	do m=2,4
800	39528000	w_left = wt_t(m) ; w_right = wt_b(m)
801	39528000	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)))
802
803		! Salinity and temperature points are linearly interpolated in
804		! the horizontal. The subscript (1) refers to the top value in
805		! the vertical profile while subscript (5) refers to the bottom
806		! value in the vertical profile.
807	39528000	pos = i15+(m-2)5
808	39528000	T15(pos+1) = (w_leftTtl) + (w_rightTtr)
809	39528000	T15(pos+5) = (w_leftTbl) + (w_rightTbr)
810
811	39528000	S15(pos+1) = (w_leftStl) + (w_rightStr)
812	39528000	S15(pos+5) = (w_leftSbl) + (w_rightSbr)
813
814	39528000	p15(pos+1) = -GxRho * ((w_left(e(i,j,K)-z0pres(i,j))) + (w_right(e(i,j+1,K)-z0pres(i,j+1))))
815
816		! Pressure
817	197640000	do n=2,5
818	197640000	p15(pos+n) = p15(pos+n-1) + GxRho0.25dz_y(m,i)
819		enddo
820
821		! Salinity and temperature (linear interpolation in the vertical)
822	158112000	do n=2,4
823	118584000	S15(pos+n) = wt_t(n) * S15(pos+1) + wt_b(n) * S15(pos+5)
824	158112000	T15(pos+n) = wt_t(n) * T15(pos+1) + wt_b(n) * T15(pos+5)
825		enddo
826	39528000	if (use_varT) T215(pos+1:pos+5) = (w_lefttv%varT(i,j,k)) + (w_righttv%varT(i,j+1,k))
827	39528000	if (use_covarTS) TS15(pos+1:pos+5) = (w_lefttv%covarTS(i,j,k)) + (w_righttv%covarTS(i,j+1,k))
828	52704000	if (use_varS) S215(pos+1:pos+5) = (w_lefttv%varS(i,j,k)) + (w_righttv%varS(i,j+1,k))
829		enddo
830		enddo
831
832	109800	if (use_stanley_eos) then
833		call calculate_density(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
834		T215(15HI%isc+1:), TS15(15HI%isc+1:), S215(15*HI%isc+1:), &
835	0	r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref)
836		else
837	109800	if (use_rho_ref) then
838		call calculate_density(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
839	109800	r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref)
840		else
841		call calculate_density(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
842	0	r15(15*HI%isc+1:), EOS, EOSdom_h15)
843		endif
844		endif
845
846	13287600	do i=HI%isc,HI%iec
847	13176000	intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
848
849		! Use Boole's rule to estimate the pressure anomaly change.
850	13176000	if (use_rho_ref) then
851	52704000	do m = 2,4
852	39528000	pos = i15+(m-2)5
853		intz(m) = (G_edz_y(m,i)( C1_90(7.0(r15(pos+1)+r15(pos+5)) + &
854		32.0*(r15(pos+2)+r15(pos+4)) + &
855	52704000	12.0*r15(pos+3)) ))
856		enddo
857		else
858	0	do m = 2,4
859	0	pos = i15+(m-2)5
860		intz(m) = (G_edz_y(m,i)( C1_90(7.0(r15(pos+1)+r15(pos+5)) + &
861		32.0*(r15(pos+2)+r15(pos+4)) + &
862	0	12.0*r15(pos+3)) - rho_ref ))
863		enddo
864		endif
865		! Use Boole's rule to integrate the values.
866		inty_dpa(i,J) = C1_90(7.0(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
867	13285800	12.0*intz(3))
868		enddo
869		enddo ; endif
870
871	7200	end subroutine int_density_dz_generic_plm
872
873
874		!> Compute pressure gradient force integrals for layer "k" and the case where T and S
875		!! are parabolic profiles
876	0	subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
877	0	rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, use_stanley_eos, &
878	0	dpa, intz_dpa, intx_dpa, inty_dpa, MassWghtInterp, Z_0p, &
879		MassWghtInterpVanOnly, h_nv)
880		integer, intent(in) :: k !< Layer index to calculate integrals for
881		type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays
882		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
883		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
884		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
885		intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC]
886		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
887		intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC]
888		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
889		intent(in) :: S_t !< Salinity at the cell top [S ~> ppt]
890		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
891		intent(in) :: S_b !< Salinity at the cell bottom [S ~> ppt]
892		real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)+1), &
893		intent(in) :: e !< Height of interfaces [Z ~> m]
894		real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is
895		!! subtracted out to reduce the magnitude of each of the integrals.
896		real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used to calculate
897		!! the pressure (as p~=-zrho_0G_e) used in the equation of state.
898		real, intent(in) :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
899		real, intent(in) :: dz_subroundoff !< A minuscule thickness change [Z ~> m]
900		real, dimension(SZI_(HI),SZJ_(HI)), &
901		intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]
902		type(EOS_type), intent(in) :: EOS !< Equation of state structure
903		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
904		logical, intent(in) :: use_stanley_eos !< If true, turn on Stanley SGS T variance parameterization
905		real, dimension(SZI_(HI),SZJ_(HI)), &
906		intent(inout) :: dpa !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]
907		real, dimension(SZI_(HI),SZJ_(HI)), &
908		optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer of
909		!! the pressure anomaly relative to the anomaly at the
910		!! top of the layer [R L2 Z T-2 ~> Pa m]
911		real, dimension(SZIB_(HI),SZJ_(HI)), &
912		optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
913		!! pressure anomaly at the top and bottom of the layer
914		!! divided by the x grid spacing [R L2 T-2 ~> Pa]
915		real, dimension(SZI_(HI),SZJB_(HI)), &
916		optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
917		!! pressure anomaly at the top and bottom of the layer
918		!! divided by the y grid spacing [R L2 T-2 ~> Pa]
919		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
920		!! mass weighting to interpolate T/S in integrals
921		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
922		!! of T/S unless one side smaller than h_nv (i.e. vanished)
923		real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]
924
925		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
926		optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
927
928		! This subroutine calculates (by numerical quadrature) integrals of
929		! pressure anomalies across layers, which are required for calculating the
930		! finite-volume form pressure accelerations in a Boussinesq model. The one
931		! potentially dodgy assumption here is that rho_0 is used both in the denominator
932		! of the accelerations, and in the pressure used to calculated density (the
933		! latter being -zrho_0G_e). These two uses could be separated if need be.
934		!
935		! It is assumed that the salinity and temperature profiles are parabolic in the
936		! vertical. The top and bottom values within each layer are provided and
937		! a parabolic interpolation is used to compute intermediate values.
938
939		! Local variables
940	0	real :: T5((5HI%iscB+1):(5(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC]
941	0	real :: S5((5HI%iscB+1):(5(HI%iecB+2))) ! Salinities along a line of subgrid locations [S ~> ppt]
942	0	real :: T25((5HI%iscB+1):(5(HI%iecB+2))) ! SGS temperature variance along a line of subgrid
943		! locations [C2 ~> degC2]
944	0	real :: TS5((5HI%iscB+1):(5(HI%iecB+2))) ! SGS temp-salt covariance along a line of subgrid
945		! locations [C S ~> degC ppt]
946	0	real :: S25((5HI%iscB+1):(5(HI%iecB+2))) ! SGS salinity variance along a line of subgrid locations [S2 ~> ppt2]
947	0	real :: p5((5HI%iscB+1):(5(HI%iecB+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]
948	0	real :: r5((5HI%iscB+1):(5(HI%iecB+2))) ! Densities anomalies along a line of subgrid
949		! locations [R ~> kg m-3]
950	0	real :: T15((15HI%iscB+1):(15(HI%iecB+1))) ! Temperatures at an array of subgrid locations [C ~> degC]
951	0	real :: S15((15HI%iscB+1):(15(HI%iecB+1))) ! Salinities at an array of subgrid locations [S ~> ppt]
952	0	real :: T215((15HI%iscB+1):(15(HI%iecB+1))) ! SGS temperature variance along a line of subgrid
953		! locations [C2 ~> degC2]
954	0	real :: TS15((15HI%iscB+1):(15(HI%iecB+1))) ! SGS temp-salt covariance along a line of subgrid
955		! locations [C S ~> degC ppt]
956	0	real :: S215((15HI%iscB+1):(15(HI%iecB+1))) ! SGS salinity variance along a line of subgrid
957		! locations [S2 ~> ppt2]
958	0	real :: p15((15HI%iscB+1):(15(HI%iecB+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]
959	0	real :: r15((15HI%iscB+1):(15(HI%iecB+1))) ! Densities at an array of subgrid locations [R ~> kg m-3]
960		real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim]
961		real :: rho_anom ! The integrated density anomaly [R ~> kg m-3]
962		real :: w_left, w_right ! Left and right weights [nondim]
963		real :: intz(5) ! The gravitational acceleration times the integrals of density
964		! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]
965		real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
966		real :: GxRho ! The gravitational acceleration times density [R L2 Z-1 T-2 ~> kg m-2 s-2]
967		real :: dz ! Layer thicknesses at tracer points [Z ~> m]
968	0	real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]
969	0	real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]
970		real :: massWeightToggle ! A non-dimensional toggle factor for near-bottom mass weighting (0 or 1) [nondim]
971		real :: TopWeightToggle ! A non-dimensional toggle factor for near-surface mass weighting (0 or 1) [nondim]
972		real :: massWeightNVonlyToggle ! A non-dimensional toggle factor for only using mass weighting
973		! if at least one side vanished (0 or 1) [nondim]
974		real :: Ttl, Tbl, Tml, Ttr, Tbr, Tmr ! Temperatures at the velocity cell corners [C ~> degC]
975		real :: Stl, Sbl, Sml, Str, Sbr, Smr ! Salinities at the velocity cell corners [S ~> ppt]
976		real :: s6 ! PPM curvature coefficient for S [S ~> ppt]
977		real :: t6 ! PPM curvature coefficient for T [C ~> degC]
978		real :: T_top, T_mn, T_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T [C ~> degC]
979		real :: S_top, S_mn, S_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S [S ~> ppt]
980	0	real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
981		real :: hWght ! A topographically limited thickness weight [Z ~> m]
982		real :: hWghtTop ! A surface displacement limited thickness weight [Z ~> m]
983		real :: hL, hR ! Thicknesses to the left and right [Z ~> m]
984		real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
985		real :: h_nonvanished ! nonvanished height [Z ~> m]
986		integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state
987		integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state
988		integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state
989		integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos
990		logical :: use_PPM ! If false, assume zero curvature in reconstruction, i.e. PLM
991		logical :: use_varT, use_varS, use_covarTS ! Logicals for SGS variances fields
992
993	0	Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB
994
995	0	GxRho = G_e * rho_0
996	0	if (present(Z_0p)) then
997	0	do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
998	0	z0pres(i,j) = Z_0p(i,j)
999		enddo ; enddo
1000		else
1001	0	z0pres(:,:) = 0.0
1002		endif
1003	0	massWeightToggle = 0. ; TopWeightToggle = 0.
1004	0	if (present(MassWghtInterp)) then
1005	0	if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1.
1006	0	if (BTEST(MassWghtInterp, 1)) TopWeightToggle = 1.
1007		endif
1008	0	massWeightNVonlyToggle = 1.
1009	0	if (present(MassWghtInterpVanOnly)) then
1010	0	if (MassWghtInterpVanOnly) massWeightNVonlyToggle = 0.
1011		endif
1012	0	h_nonvanished = 0.
1013	0	if (present(h_nv)) then
1014	0	h_nonvanished = h_nv
1015		endif
1016
1017		! In event PPM calculation is bypassed with use_PPM=False
1018	0	s6 = 0.
1019	0	t6 = 0.
1020	0	use_PPM = .true. ! This is a place-holder to allow later re-use of this function
1021
1022	0	use_varT = .false. !ensure initialized
1023	0	use_covarTS = .false.
1024	0	use_varS = .false.
1025	0	if (use_stanley_eos) then
1026	0	use_varT = associated(tv%varT)
1027	0	use_covarTS = associated(tv%covarTS)
1028	0	use_varS = associated(tv%varS)
1029		endif
1030
1031	0	T25(:) = 0.
1032	0	TS5(:) = 0.
1033	0	S25(:) = 0.
1034	0	T215(:) = 0.
1035	0	TS15(:) = 0.
1036	0	S215(:) = 0.
1037
1038	0	do n = 1, 5
1039	0	wt_t(n) = 0.25 * real(5-n)
1040	0	wt_b(n) = 1.0 - wt_t(n)
1041		enddo
1042
1043		! Set the loop ranges for equation of state calculations at various points.
1044	0	EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2)
1045	0	EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1)
1046	0	EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1)
1047
1048		! 1. Compute vertical integrals
1049	0	do j=Jsq,Jeq+1
1050	0	do i=Isq,Ieq+1
1051	0	if (use_PPM) then
1052		! Curvature coefficient of the parabolas
1053	0	s6 = 3.0 * ( 2.0*tv%S(i,j,k) - ( S_t(i,j,k) + S_b(i,j,k) ) )
1054	0	t6 = 3.0 * ( 2.0*tv%T(i,j,k) - ( T_t(i,j,k) + T_b(i,j,k) ) )
1055		endif
1056	0	dz = e(i,j,K) - e(i,j,K+1)
1057	0	do n=1,5
1058	0	p5(I5+n) = -GxRho((e(i,j,K) - z0pres(i,j)) - 0.25real(n-1)dz)
1059		! Salinity and temperature points are reconstructed with PPM
1060	0	S5(I5+n) = wt_t(n) S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) )
1061	0	T5(I5+n) = wt_t(n) T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) )
1062		enddo
1063	0	if (use_stanley_eos) then
1064	0	if (use_varT) T25(I5+1:I5+5) = tv%varT(i,j,k)
1065	0	if (use_covarTS) TS5(I5+1:I5+5) = tv%covarTS(i,j,k)
1066	0	if (use_varS) S25(I5+1:I5+5) = tv%varS(i,j,k)
1067		endif
1068		enddo
1069
1070	0	if (use_stanley_eos) then
1071	0	call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, EOSdom_h5, rho_ref=rho_ref)
1072		else
1073	0	call calculate_density(T5, S5, p5, r5, EOS, EOSdom_h5, rho_ref=rho_ref)
1074		endif
1075
1076	0	do i=Isq,Ieq+1
1077	0	dz = e(i,j,K) - e(i,j,K+1)
1078		! Use Boole's rule to estimate the pressure anomaly change.
1079	0	rho_anom = C1_90(7.0(r5(i5+1)+r5(i5+5)) + 32.0(r5(i5+2)+r5(i5+4)) + 12.0r5(i*5+3))
1080	0	dpa(i,j) = G_edzrho_anom
1081	0	if (present(intz_dpa)) then
1082		! Use a Boole's-rule-like fifth-order accurate estimate of
1083		! the double integral of the pressure anomaly.
1084		intz_dpa(i,j) = 0.5G_edz*2 &
1085	0	(rho_anom - C1_90(16.0(r5(i5+4)-r5(i5+2)) + 7.0(r5(i5+5)-r5(i*5+1))) )
1086		endif
1087		enddo ! end loop on i
1088		enddo ! end loop on j
1089
1090		! 2. Compute horizontal integrals in the x direction
1091	0	if (present(intx_dpa)) then ; do j=HI%jsc,HI%jec
1092	0	do I=Isq,Ieq
1093		! Corner values of T and S
1094		! hWght is the distance measure by which the cell is violation of
1095		! hydrostatic consistency. For large hWght we bias the interpolation
1096		! of T,S along the top and bottom integrals, almost like thickness
1097		! weighting.
1098		! Note: To work in terrain following coordinates we could offset
1099		! this distance by the layer thickness to replicate other models.
1100		hWght = massWeightToggle * &
1101	0	max(0., -bathyT(i,j)-e(i+1,j,K), -bathyT(i+1,j)-e(i,j,K))
1102		hWghtTop = TopWeightToggle * &
1103	0	max(0., e(i+1,j,K+1)-e(i,j,1), e(i,j,K+1)-e(i+1,j,1))
1104	0	hWght = max(hWght, hWghtTop)
1105		! If both sides are nonvanished, then set it back to zero.
1106	0	if (((e(i,j,K) - e(i,j,K+1)) > h_nonvanished) .and. ((e(i+1,j,K) - e(i+1,j,K+1)) > h_nonvanished)) then
1107	0	hWght = massWeightNVonlyToggle * hWght
1108		endif
1109	0	if (hWght > 0.) then
1110	0	hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff
1111	0	hR = (e(i+1,j,K) - e(i+1,j,K+1)) + dz_subroundoff
1112	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
1113	0	iDenom = 1./( hWght(hR + hL) + hLhR )
1114	0	Ttl = ( (hWghthR)T_t(i+1,j,k) + (hWghthL + hRhL)T_t(i,j,k) ) iDenom
1115	0	Tbl = ( (hWghthR)T_b(i+1,j,k) + (hWghthL + hRhL)T_b(i,j,k) ) iDenom
1116	0	Tml = ( (hWghthR)tv%T(i+1,j,k)+ (hWghthL + hRhL)tv%T(i,j,k) ) iDenom
1117	0	Ttr = ( (hWghthL)T_t(i,j,k) + (hWghthR + hRhL)T_t(i+1,j,k) ) iDenom
1118	0	Tbr = ( (hWghthL)T_b(i,j,k) + (hWghthR + hRhL)T_b(i+1,j,k) ) iDenom
1119	0	Tmr = ( (hWghthL)tv%T(i,j,k) + (hWghthR + hRhL)tv%T(i+1,j,k) ) iDenom
1120	0	Stl = ( (hWghthR)S_t(i+1,j,k) + (hWghthL + hRhL)S_t(i,j,k) ) iDenom
1121	0	Sbl = ( (hWghthR)S_b(i+1,j,k) + (hWghthL + hRhL)S_b(i,j,k) ) iDenom
1122	0	Sml = ( (hWghthR)tv%S(i+1,j,k) + (hWghthL + hRhL)tv%S(i,j,k) ) iDenom
1123	0	Str = ( (hWghthL)S_t(i,j,k) + (hWghthR + hRhL)S_t(i+1,j,k) ) iDenom
1124	0	Sbr = ( (hWghthL)S_b(i,j,k) + (hWghthR + hRhL)S_b(i+1,j,k) ) iDenom
1125	0	Smr = ( (hWghthL)tv%S(i,j,k) + (hWghthR + hRhL)tv%S(i+1,j,k) ) iDenom
1126		else
1127	0	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)
1128	0	Tml = tv%T(i,j,k) ; Tmr = tv%T(i+1,j,k)
1129	0	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)
1130	0	Sml = tv%S(i,j,k) ; Smr = tv%S(i+1,j,k)
1131		endif
1132
1133	0	do m=2,4
1134	0	w_left = wt_t(m) ; w_right = wt_b(m)
1135
1136		! Salinity and temperature points are linearly interpolated in
1137		! the horizontal. The subscript (1) refers to the top value in
1138		! the vertical profile while subscript (5) refers to the bottom
1139		! value in the vertical profile.
1140	0	T_top = (w_leftTtl) + (w_rightTtr)
1141	0	T_mn = (w_leftTml) + (w_rightTmr)
1142	0	T_bot = (w_leftTbl) + (w_rightTbr)
1143
1144	0	S_top = (w_leftStl) + (w_rightStr)
1145	0	S_mn = (w_leftSml) + (w_rightSmr)
1146	0	S_bot = (w_leftSbl) + (w_rightSbr)
1147
1148		! Pressure
1149	0	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)))
1150
1151	0	pos = i15+(m-2)5
1152	0	p15(pos+1) = -GxRho * ((w_left(e(i,j,K)-z0pres(i,j))) + (w_right(e(i+1,j,K)-z0pres(i+1,j))))
1153	0	do n=2,5
1154	0	p15(pos+n) = p15(pos+n-1) + GxRho0.25dz_x(m,i)
1155		enddo
1156
1157		! Parabolic reconstructions in the vertical for T and S
1158	0	if (use_PPM) then
1159		! Coefficients of the parabolas
1160	0	s6 = 3.0 * ( 2.0*S_mn - ( S_top + S_bot ) )
1161	0	t6 = 3.0 * ( 2.0*T_mn - ( T_top + T_bot ) )
1162		endif
1163	0	do n=1,5
1164	0	S15(pos+n) = wt_t(n) * S_top + wt_b(n) * ( S_bot + s6 * wt_t(n) )
1165	0	T15(pos+n) = wt_t(n) * T_top + wt_b(n) * ( T_bot + t6 * wt_t(n) )
1166		enddo
1167	0	if (use_stanley_eos) then
1168	0	if (use_varT) T215(pos+1:pos+5) = (w_lefttv%varT(i,j,k)) + (w_righttv%varT(i+1,j,k))
1169	0	if (use_covarTS) TS15(pos+1:pos+5) = (w_lefttv%covarTS(i,j,k)) + (w_righttv%covarTS(i+1,j,k))
1170	0	if (use_varS) S215(pos+1:pos+5) = (w_lefttv%varS(i,j,k)) + (w_righttv%varS(i+1,j,k))
1171		endif
1172	0	if (use_stanley_eos) then
1173	0	call calculate_density(T5, S5, p5, T25, TS5, S25, r5, EOS, rho_ref=rho_ref)
1174		else
1175	0	call calculate_density(T5, S5, p5, r5, EOS, rho_ref=rho_ref)
1176		endif
1177		enddo
1178		enddo
1179
1180	0	if (use_stanley_eos) then
1181	0	call calculate_density(T15, S15, p15, T215, TS15, S215, r15, EOS, EOSdom_q15, rho_ref=rho_ref)
1182		else
1183	0	call calculate_density(T15, S15, p15, r15, EOS, EOSdom_q15, rho_ref=rho_ref)
1184		endif
1185
1186	0	do I=Isq,Ieq
1187	0	do m=2,4
1188	0	pos = i15+(m-2)5
1189		! Use Boole's rule to estimate the pressure anomaly change.
1190		intz(m) = (G_edz_x(m,i)(C1_90( 7.0(r15(pos+1)+r15(pos+5)) + &
1191		32.0*(r15(pos+2)+r15(pos+4)) + &
1192	0	12.0*r15(pos+3)) ))
1193		enddo ! m
1194	0	intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
1195
1196		! Use Boole's rule to integrate the bottom pressure anomaly values in x.
1197	0	intx_dpa(I,j) = C1_90(7.0(intz(1)+intz(5)) + 32.0(intz(2)+intz(4)) + 12.0intz(3))
1198
1199		enddo
1200		enddo ; endif
1201
1202		! 3. Compute horizontal integrals in the y direction
1203	0	if (present(inty_dpa)) then ; do J=Jsq,Jeq
1204	0	do i=HI%isc,HI%iec
1205		! Corner values of T and S
1206		! hWght is the distance measure by which the cell is violation of
1207		! hydrostatic consistency. For large hWght we bias the interpolation
1208		! of T,S along the top and bottom integrals, almost like thickness
1209		! weighting.
1210		! Note: To work in terrain following coordinates we could offset
1211		! this distance by the layer thickness to replicate other models.
1212		hWght = massWeightToggle * &
1213	0	max(0., -bathyT(i,j)-e(i,j+1,K), -bathyT(i,j+1)-e(i,j,K))
1214		hWghtTop = TopWeightToggle * &
1215	0	max(0., e(i,j+1,K+1)-e(i,j,1), e(i,j,K+1)-e(i,j+1,1))
1216	0	hWght = max(hWght, hWghtTop)
1217		! If both sides are nonvanished, then set it back to zero.
1218	0	if (((e(i,j,K) - e(i,j,K+1)) > h_nonvanished) .and. ((e(i,j+1,K) - e(i,j+1,K+1)) > h_nonvanished)) then
1219	0	hWght = massWeightNVonlyToggle * hWght
1220		endif
1221	0	if (hWght > 0.) then
1222	0	hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff
1223	0	hR = (e(i,j+1,K) - e(i,j+1,K+1)) + dz_subroundoff
1224	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
1225	0	iDenom = 1./( hWght(hR + hL) + hLhR )
1226	0	Ttl = ( (hWghthR)T_t(i,j+1,k) + (hWghthL + hRhL)T_t(i,j,k) ) iDenom
1227	0	Tbl = ( (hWghthR)T_b(i,j+1,k) + (hWghthL + hRhL)T_b(i,j,k) ) iDenom
1228	0	Tml = ( (hWghthR)tv%T(i,j+1,k)+ (hWghthL + hRhL)tv%T(i,j,k) ) iDenom
1229	0	Ttr = ( (hWghthL)T_t(i,j,k) + (hWghthR + hRhL)T_t(i,j+1,k) ) iDenom
1230	0	Tbr = ( (hWghthL)T_b(i,j,k) + (hWghthR + hRhL)T_b(i,j+1,k) ) iDenom
1231	0	Tmr = ( (hWghthL)tv%T(i,j,k) + (hWghthR + hRhL)tv%T(i,j+1,k) ) iDenom
1232	0	Stl = ( (hWghthR)S_t(i,j+1,k) + (hWghthL + hRhL)S_t(i,j,k) ) iDenom
1233	0	Sbl = ( (hWghthR)S_b(i,j+1,k) + (hWghthL + hRhL)S_b(i,j,k) ) iDenom
1234	0	Sml = ( (hWghthR)tv%S(i,j+1,k)+ (hWghthL + hRhL)tv%S(i,j,k) ) iDenom
1235	0	Str = ( (hWghthL)S_t(i,j,k) + (hWghthR + hRhL)S_t(i,j+1,k) ) iDenom
1236	0	Sbr = ( (hWghthL)S_b(i,j,k) + (hWghthR + hRhL)S_b(i,j+1,k) ) iDenom
1237	0	Smr = ( (hWghthL)tv%S(i,j,k) + (hWghthR + hRhL)tv%S(i,j+1,k) ) iDenom
1238		else
1239	0	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)
1240	0	Tml = tv%T(i,j,k) ; Tmr = tv%T(i,j+1,k)
1241	0	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)
1242	0	Sml = tv%S(i,j,k) ; Smr = tv%S(i,j+1,k)
1243		endif
1244
1245	0	do m=2,4
1246	0	w_left = wt_t(m) ; w_right = wt_b(m)
1247
1248		! Salinity and temperature points are linearly interpolated in
1249		! the horizontal. The subscript (1) refers to the top value in
1250		! the vertical profile while subscript (5) refers to the bottom
1251		! value in the vertical profile.
1252	0	T_top = (w_leftTtl) + (w_rightTtr)
1253	0	T_mn = (w_leftTml) + (w_rightTmr)
1254	0	T_bot = (w_leftTbl) + (w_rightTbr)
1255
1256	0	S_top = (w_leftStl) + (w_rightStr)
1257	0	S_mn = (w_leftSml) + (w_rightSmr)
1258	0	S_bot = (w_leftSbl) + (w_rightSbr)
1259
1260		! Pressure
1261	0	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)))
1262
1263	0	pos = i15+(m-2)5
1264	0	p15(pos+1) = -GxRho * ((w_left(e(i,j,K)-z0pres(i,j))) + (w_right(e(i,j+1,K)-z0pres(i,j+1))))
1265	0	do n=2,5
1266	0	p15(pos+n) = p15(pos+n-1) + GxRho0.25dz_y(m,i)
1267		enddo
1268
1269		! Parabolic reconstructions in the vertical for T and S
1270	0	if (use_PPM) then
1271		! Coefficients of the parabolas
1272	0	s6 = 3.0 * ( 2.0*S_mn - ( S_top + S_bot ) )
1273	0	t6 = 3.0 * ( 2.0*T_mn - ( T_top + T_bot ) )
1274		endif
1275	0	do n=1,5
1276	0	S15(pos+n) = wt_t(n) * S_top + wt_b(n) * ( S_bot + s6 * wt_t(n) )
1277	0	T15(pos+n) = wt_t(n) * T_top + wt_b(n) * ( T_bot + t6 * wt_t(n) )
1278		enddo
1279
1280	0	if (use_stanley_eos) then
1281	0	if (use_varT) T215(pos+1:pos+5) = (w_lefttv%varT(i,j,k)) + (w_righttv%varT(i,j+1,k))
1282	0	if (use_covarTS) TS15(pos+1:pos+5) = (w_lefttv%covarTS(i,j,k)) + (w_righttv%covarTS(i,j+1,k))
1283	0	if (use_varS) S215(pos+1:pos+5) = (w_lefttv%varS(i,j,k)) + (w_righttv%varS(i,j+1,k))
1284		endif
1285		enddo
1286		enddo
1287
1288	0	if (use_stanley_eos) then
1289		call calculate_density(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
1290		T215(15HI%isc+1:), TS15(15HI%isc+1:), S215(15*HI%isc+1:), &
1291	0	r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref)
1292		else
1293		call calculate_density(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
1294	0	r15(15*HI%isc+1:), EOS, EOSdom_h15, rho_ref=rho_ref)
1295		endif
1296
1297	0	do i=HI%isc,HI%iec
1298	0	do m=2,4
1299		! Use Boole's rule to estimate the pressure anomaly change.
1300	0	pos = i15+(m-2)5
1301		intz(m) = (G_edz_y(m,i)(C1_90( 7.0(r15(pos+1)+r15(pos+5)) + &
1302		32.0*(r15(pos+2)+r15(pos+4)) + &
1303	0	12.0*r15(pos+3)) ))
1304		enddo ! m
1305	0	intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
1306
1307		! Use Boole's rule to integrate the bottom pressure anomaly values in y.
1308	0	inty_dpa(i,J) = C1_90(7.0(intz(1)+intz(5)) + 32.0(intz(2)+intz(4)) + 12.0intz(3))
1309		enddo
1310		enddo ; endif
1311
1312	0	end subroutine int_density_dz_generic_ppm
1313
1314		!> Calls the appropriate subroutine to calculate analytical and nearly-analytical
1315		!! integrals in pressure across layers of geopotential anomalies, which are
1316		!! required for calculating the finite-volume form pressure accelerations in a
1317		!! non-Boussinesq model. There are essentially no free assumptions, apart from the
1318		!! use of Boole's rule to do the horizontal integrals, and from a truncation in the
1319		!! series for log(1-eps/1+eps) that assumes that \|eps\| < 0.34.
1320	0	subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, &
1321	0	dza, intp_dza, intx_dza, inty_dza, halo_size, &
1322	0	bathyP, P_surf, dP_tiny, MassWghtInterp, &
1323		MassWghtInterpVanOnly, p_nv)
1324		type(hor_index_type), intent(in) :: HI !< The horizontal index structure
1325		real, dimension(SZI_(HI),SZJ_(HI)), &
1326		intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
1327		real, dimension(SZI_(HI),SZJ_(HI)), &
1328		intent(in) :: S !< Salinity [S ~> ppt]
1329		real, dimension(SZI_(HI),SZJ_(HI)), &
1330		intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
1331		real, dimension(SZI_(HI),SZJ_(HI)), &
1332		intent(in) :: p_b !< Pressure at the bottom of the layer [R L2 T-2 ~> Pa]
1333		real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
1334		!! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
1335		!! The calculation is mathematically identical with different values of
1336		!! alpha_ref, but this reduces the effects of roundoff.
1337		type(EOS_type), intent(in) :: EOS !< Equation of state structure
1338		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1339		real, dimension(SZI_(HI),SZJ_(HI)), &
1340		intent(inout) :: dza !< The change in the geopotential anomaly across
1341		!! the layer [L2 T-2 ~> m2 s-2]
1342		real, dimension(SZI_(HI),SZJ_(HI)), &
1343		optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the
1344		!! geopotential anomaly relative to the anomaly at the bottom of the
1345		!! layer [R L4 T-4 ~> Pa m2 s-2]
1346		real, dimension(SZIB_(HI),SZJ_(HI)), &
1347		optional, intent(inout) :: intx_dza !< The integral in x of the difference between the
1348		!! geopotential anomaly at the top and bottom of the layer divided by
1349		!! the x grid spacing [L2 T-2 ~> m2 s-2]
1350		real, dimension(SZI_(HI),SZJB_(HI)), &
1351		optional, intent(inout) :: inty_dza !< The integral in y of the difference between the
1352		!! geopotential anomaly at the top and bottom of the layer divided by
1353		!! the y grid spacing [L2 T-2 ~> m2 s-2]
1354		integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
1355		real, dimension(SZI_(HI),SZJ_(HI)), &
1356		optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
1357		real, dimension(SZI_(HI),SZJ_(HI)), &
1358		optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
1359		real, optional, intent(in) :: dP_tiny !< A minuscule pressure change with
1360		!! the same units as p_t [R L2 T-2 ~> Pa]
1361		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
1362		!! mass weighting to interpolate T/S in integrals
1363		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
1364		!! of T/S unless one side smaller than h_nv (i.e. vanished)
1365		real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]
1366
1367
1368	0	if (EOS_quadrature(EOS)) then
1369		call int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, &
1370		dza, intp_dza, intx_dza, inty_dza, halo_size, &
1371		bathyP, P_surf, dP_tiny, MassWghtInterp, &
1372	0	MassWghtInterpVanOnly, p_nv)
1373		else
1374		call analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, &
1375		dza, intp_dza, intx_dza, inty_dza, halo_size, &
1376	0	bathyP, P_surf, dP_tiny, MassWghtInterp)
1377		endif
1378
1379	0	end subroutine int_specific_vol_dp
1380
1381
1382		!> This subroutine calculates integrals of specific volume anomalies in
1383		!! pressure across layers, which are required for calculating the finite-volume
1384		!! form pressure accelerations in a non-Boussinesq model. There are essentially
1385		!! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.
1386	0	subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, dza, &
1387	0	intp_dza, intx_dza, inty_dza, halo_size, &
1388	0	bathyP, P_surf, dP_neglect, MassWghtInterp, &
1389		MassWghtInterpVanOnly, p_nv)
1390		type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
1391		real, dimension(SZI_(HI),SZJ_(HI)), &
1392		intent(in) :: T !< Potential temperature of the layer [C ~> degC]
1393		real, dimension(SZI_(HI),SZJ_(HI)), &
1394		intent(in) :: S !< Salinity of the layer [S ~> ppt]
1395		real, dimension(SZI_(HI),SZJ_(HI)), &
1396		intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa]
1397		real, dimension(SZI_(HI),SZJ_(HI)), &
1398		intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa]
1399		real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
1400		!! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
1401		!! The calculation is mathematically identical with different values of
1402		!! alpha_ref, but alpha_ref alters the effects of roundoff, and
1403		!! answers do change.
1404		type(EOS_type), intent(in) :: EOS !< Equation of state structure
1405		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1406		real, dimension(SZI_(HI),SZJ_(HI)), &
1407		intent(inout) :: dza !< The change in the geopotential anomaly
1408		!! across the layer [L2 T-2 ~> m2 s-2]
1409		real, dimension(SZI_(HI),SZJ_(HI)), &
1410		optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
1411		!! the geopotential anomaly relative to the anomaly at the bottom of the
1412		!! layer [R L4 T-4 ~> Pa m2 s-2]
1413		real, dimension(SZIB_(HI),SZJ_(HI)), &
1414		optional, intent(inout) :: intx_dza !< The integral in x of the difference between
1415		!! the geopotential anomaly at the top and bottom of the layer divided
1416		!! by the x grid spacing [L2 T-2 ~> m2 s-2]
1417		real, dimension(SZI_(HI),SZJB_(HI)), &
1418		optional, intent(inout) :: inty_dza !< The integral in y of the difference between
1419		!! the geopotential anomaly at the top and bottom of the layer divided
1420		!! by the y grid spacing [L2 T-2 ~> m2 s-2]
1421		integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
1422		real, dimension(SZI_(HI),SZJ_(HI)), &
1423		optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
1424		real, dimension(SZI_(HI),SZJ_(HI)), &
1425		optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
1426		real, optional, intent(in) :: dP_neglect !< A minuscule pressure change with
1427		!! the same units as p_t [R L2 T-2 ~> Pa]
1428		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
1429		!! mass weighting to interpolate T/S in integrals
1430		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
1431		!! of T/S unless one side smaller than h_nv (i.e. vanished)
1432		real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]
1433
1434		! This subroutine calculates analytical and nearly-analytical integrals in
1435		! pressure across layers of geopotential anomalies, which are required for
1436		! calculating the finite-volume form pressure accelerations in a non-Boussinesq
1437		! model. There are essentially no free assumptions, apart from the use of
1438		! Boole's rule to do the horizontal integrals, and from a truncation in the
1439		! series for log(1-eps/1+eps) that assumes that \|eps\| < 0.34.
1440
1441		! Local variables
1442	0	real :: T5((5HI%isd+1):(5(HI%ied+2))) ! Temperatures along a line of subgrid locations [C ~> degC]
1443	0	real :: S5((5HI%isd+1):(5(HI%ied+2))) ! Salinities along a line of subgrid locations [S ~> ppt]
1444	0	real :: p5((5HI%isd+1):(5(HI%ied+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]
1445	0	real :: a5((5HI%isd+1):(5(HI%ied+2))) ! Specific volumes anomalies along a line of subgrid
1446		! locations [R-1 ~> m3 kg-1]
1447	0	real :: T15((15HI%isd+1):(15(HI%ied+1))) ! Temperatures at an array of subgrid locations [C ~> degC]
1448	0	real :: S15((15HI%isd+1):(15(HI%ied+1))) ! Salinities at an array of subgrid locations [S ~> ppt]
1449	0	real :: p15((15HI%isd+1):(15(HI%ied+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]
1450	0	real :: a15((15HI%isd+1):(15(HI%ied+1))) ! Specific volumes at an array of subgrid locations [R ~> kg m-3]
1451		real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]
1452		real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]
1453	0	real :: dp_x(5,SZIB_(HI)) ! The pressure change through a layer along an x-line of subgrid locations [Z ~> m]
1454	0	real :: dp_y(5,SZI_(HI)) ! The pressure change through a layer along a y-line of subgrid locations [Z ~> m]
1455		real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
1456		real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
1457		real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
1458		real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim]
1459		real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]
1460		real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]
1461		real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]
1462		real :: intp(5) ! The integrals of specific volume with pressure at the
1463		! 5 sub-column locations [L2 T-2 ~> m2 s-2]
1464		logical :: do_massWeight ! Indicates whether to do mass weighting.
1465		logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
1466		logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting
1467		real :: massWeightNVonlyToggle ! A non-dimensional toggle factor for only using mass weighting
1468		! if at least one side vanished (0 or 1) [nondim]
1469		real :: p_nonvanished ! nonvanished pressure [R L2 T-2 ~> Pa]
1470		real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
1471		integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state
1472		integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state
1473		integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state
1474		integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, n, pos, halo
1475
1476	0	Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB
1477	0	halo = 0 ; if (present(halo_size)) halo = MAX(halo_size,0)
1478	0	ish = HI%isc-halo ; ieh = HI%iec+halo ; jsh = HI%jsc-halo ; jeh = HI%jec+halo
1479	0	if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh) ; endif
1480	0	if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh) ; endif
1481
1482	0	do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false.
1483	0	if (present(MassWghtInterp)) then
1484	0	do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
1485	0	top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
1486	0	massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set
1487	0	if (do_massWeight .and. .not.present(bathyP)) call MOM_error(FATAL, &
1488	0	"int_spec_vol_dp_generic_pcm: bathyP must be present if near-bottom mass weighting is in use.")
1489	0	if (top_massWeight .and. .not.present(P_surf)) call MOM_error(FATAL, &
1490	0	"int_spec_vol_dp_generic_pcm: P_surf must be present if near-surface mass weighting is in use.")
1491	0	if ((do_massWeight .or. top_massWeight) .and. .not.present(dP_neglect)) call MOM_error(FATAL, &
1492	0	"int_spec_vol_dp_generic_pcm: dP_neglect must be present if mass weighting is in use.")
1493		endif
1494	0	massWeightNVonlyToggle = 1.
1495	0	if (present(MassWghtInterpVanOnly)) then
1496	0	if (MassWghtInterpVanOnly) massWeightNVonlyToggle = 0.
1497		endif
1498	0	p_nonvanished = 0.
1499	0	if (present(p_nv)) then
1500	0	p_nonvanished = p_nv
1501		endif
1502
1503
1504		! Set the loop ranges for equation of state calculations at various points.
1505	0	EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(ieh-ish+1)
1506	0	EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1)
1507	0	EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1)
1508
1509	0	do j=jsh,jeh
1510	0	do i=ish,ieh
1511	0	dp = p_b(i,j) - p_t(i,j)
1512	0	pos = 5*i
1513	0	do n=1,5
1514	0	T5(pos+n) = T(i,j) ; S5(pos+n) = S(i,j)
1515	0	p5(pos+n) = p_b(i,j) - 0.25real(n-1)dp
1516		enddo
1517		enddo
1518
1519		call calculate_spec_vol(T5(5ish+1:), S5(5ish+1:), p5(5ish+1:), a5(5ish+1:), EOS, &
1520	0	EOSdom_h5, spv_ref=alpha_ref)
1521
1522	0	do i=ish,ieh
1523	0	dp = p_b(i,j) - p_t(i,j)
1524		! Use Boole's rule to estimate the interface height anomaly change.
1525	0	pos = 5*i
1526	0	alpha_anom = C1_90(7.0(a5(pos+1)+a5(pos+5)) + 32.0(a5(pos+2)+a5(pos+4)) + 12.0a5(pos+3))
1527	0	dza(i,j) = dp*alpha_anom
1528		! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
1529		! the interface height anomaly.
1530	0	if (present(intp_dza)) intp_dza(i,j) = 0.5dp2 &
1531	0	(alpha_anom - C1_90(16.0(a5(pos+4)-a5(pos+2)) + 7.0*(a5(pos+5)-a5(pos+1))) )
1532		enddo
1533		enddo
1534
1535	0	if (present(intx_dza)) then ; do j=HI%jsc,HI%jec
1536	0	do I=Isq,Ieq
1537		! hWght is the distance measure by which the cell is violation of
1538		! hydrostatic consistency. For large hWght we bias the interpolation of
1539		! T & S along the top and bottom integrals, akin to thickness weighting.
1540	0	hWght = 0.0
1541	0	if (do_massWeight .and. massWeight_bug) then
1542	0	hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j))
1543	0	elseif (do_massWeight) then
1544	0	hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j))
1545		endif
1546	0	if (top_massWeight) &
1547	0	hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j))
1548		! If both sides are nonvanished, then set it back to zero.
1549	0	if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i+1,j) - p_t(i+1,j)) > p_nonvanished)) then
1550	0	hWght = massWeightNVonlyToggle * hWght
1551		endif
1552
1553	0	if (hWght > 0.) then
1554	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
1555	0	hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
1556	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
1557	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
1558	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
1559	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
1560		else
1561	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
1562		endif
1563
1564	0	do m=2,4
1565	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
1566	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
1567	0	pos = i15+(m-2)5
1568
1569		! T, S, and p are interpolated in the horizontal. The p interpolation
1570		! is linear, but for T and S it may be thickness weighted.
1571	0	p15(pos+1) = (wt_Lp_b(i,j)) + (wt_Rp_b(i+1,j))
1572	0	dp_x(m,I) = (wt_L(p_b(i,j) - p_t(i,j))) + (wt_R(p_b(i+1,j) - p_t(i+1,j)))
1573	0	T15(pos+1) = (wtT_LT(i,j)) + (wtT_RT(i+1,j))
1574	0	S15(pos+1) = (wtT_LS(i,j)) + (wtT_RS(i+1,j))
1575
1576	0	do n=2,5
1577	0	T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
1578	0	p15(pos+n) = p15(pos+n-1) - 0.25*dp_x(m,I)
1579		enddo
1580		enddo
1581		enddo
1582
1583		call calculate_spec_vol(T15(15Isq+1:), S15(15Isq+1:), p15(15*Isq+1:), &
1584	0	a15(15*Isq+1:), EOS, EOSdom_q15, spv_ref=alpha_ref)
1585
1586	0	do I=Isq,Ieq
1587	0	intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
1588		! Use Boole's rule to estimate the interface height anomaly change.
1589	0	do m=2,4
1590	0	pos = i15+(m-2)5
1591		intp(m) = (dp_x(m,I)( C1_90(7.0(a15(pos+1)+a15(pos+5)) + 32.0(a15(pos+2)+a15(pos+4)) + &
1592	0	12.0*a15(pos+3)) ))
1593		enddo
1594		! Use Boole's rule to integrate the interface height anomaly values in x.
1595		intx_dza(i,j) = C1_90(7.0(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
1596	0	12.0*intp(3))
1597		enddo
1598		enddo ; endif
1599
1600	0	if (present(inty_dza)) then ; do J=Jsq,Jeq
1601	0	do i=HI%isc,HI%iec
1602		! hWght is the distance measure by which the cell is violation of
1603		! hydrostatic consistency. For large hWght we bias the interpolation of
1604		! T & S along the top and bottom integrals, akin to thickness weighting.
1605	0	hWght = 0.0
1606	0	if (do_massWeight .and. massWeight_bug) then
1607	0	hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j))
1608	0	elseif (do_massWeight) then
1609	0	hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1))
1610		endif
1611	0	if (top_massWeight) &
1612	0	hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j))
1613		! If both sides are nonvanished, then set it back to zero.
1614	0	if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i,j+1) - p_t(i,j+1)) > p_nonvanished)) then
1615	0	hWght = massWeightNVonlyToggle * hWght
1616		endif
1617	0	if (hWght > 0.) then
1618	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
1619	0	hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
1620	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
1621	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
1622	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
1623	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
1624		else
1625	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
1626		endif
1627
1628	0	do m=2,4
1629	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
1630	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
1631	0	pos = i15+(m-2)5
1632
1633		! T, S, and p are interpolated in the horizontal. The p interpolation
1634		! is linear, but for T and S it may be thickness weighted.
1635	0	p15(pos+1) = (wt_Lp_b(i,j)) + (wt_Rp_b(i,j+1))
1636	0	dp_y(m,i) = (wt_L(p_b(i,j) - p_t(i,j))) + (wt_R(p_b(i,j+1) - p_t(i,j+1)))
1637	0	T15(pos+1) = (wtT_LT(i,j)) + (wtT_RT(i,j+1))
1638	0	S15(pos+1) = (wtT_LS(i,j)) + (wtT_RS(i,j+1))
1639	0	do n=2,5
1640	0	T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
1641	0	p15(pos+n) = p15(pos+n-1) - 0.25*dp_y(m,i)
1642		enddo
1643		enddo
1644		enddo
1645
1646		call calculate_spec_vol(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
1647	0	a15(15*HI%isc+1:), EOS, EOSdom_h15, spv_ref=alpha_ref)
1648
1649	0	do i=HI%isc,HI%iec
1650
1651	0	intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
1652		! Use Boole's rule to estimate the interface height anomaly change.
1653	0	do m=2,4
1654	0	pos = i15+(m-2)5
1655		intp(m) = (dp_y(m,i)( C1_90(7.0(a15(pos+1)+a15(pos+5)) + 32.0(a15(pos+2)+a15(pos+4)) + &
1656	0	12.0*a15(pos+3)) ))
1657		enddo
1658		! Use Boole's rule to integrate the interface height anomaly values in y.
1659		inty_dza(i,j) = C1_90(7.0(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
1660	0	12.0*intp(3))
1661		enddo
1662		enddo ; endif
1663
1664	0	end subroutine int_spec_vol_dp_generic_pcm
1665
1666		!> This subroutine calculates integrals of specific volume anomalies in
1667		!! pressure across layers, which are required for calculating the finite-volume
1668		!! form pressure accelerations in a non-Boussinesq model. There are essentially
1669		!! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals.
1670	0	subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, &
1671	0	dP_neglect, bathyP, HI, EOS, US, dza, &
1672	0	intp_dza, intx_dza, inty_dza, P_surf, MassWghtInterp, &
1673		MassWghtInterpVanOnly, p_nv)
1674		type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
1675		real, dimension(SZI_(HI),SZJ_(HI)), &
1676		intent(in) :: T_t !< Potential temperature at the top of the layer [C ~> degC]
1677		real, dimension(SZI_(HI),SZJ_(HI)), &
1678		intent(in) :: T_b !< Potential temperature at the bottom of the layer [C ~> degC]
1679		real, dimension(SZI_(HI),SZJ_(HI)), &
1680		intent(in) :: S_t !< Salinity at the top the layer [S ~> ppt]
1681		real, dimension(SZI_(HI),SZJ_(HI)), &
1682		intent(in) :: S_b !< Salinity at the bottom the layer [S ~> ppt]
1683		real, dimension(SZI_(HI),SZJ_(HI)), &
1684		intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa]
1685		real, dimension(SZI_(HI),SZJ_(HI)), &
1686		intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa]
1687		real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
1688		!! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
1689		!! The calculation is mathematically identical with different values of
1690		!! alpha_ref, but alpha_ref alters the effects of roundoff, and
1691		!! answers do change.
1692		real, intent(in) :: dP_neglect !<!< A miniscule pressure change with
1693		!! the same units as p_t [R L2 T-2 ~> Pa]
1694		real, dimension(SZI_(HI),SZJ_(HI)), &
1695		intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
1696		type(EOS_type), intent(in) :: EOS !< Equation of state structure
1697		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1698		real, dimension(SZI_(HI),SZJ_(HI)), &
1699		intent(inout) :: dza !< The change in the geopotential anomaly
1700		!! across the layer [L2 T-2 ~> m2 s-2]
1701		real, dimension(SZI_(HI),SZJ_(HI)), &
1702		optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
1703		!! the geopotential anomaly relative to the anomaly at the bottom of the
1704		!! layer [R L4 T-4 ~> Pa m2 s-2]
1705		real, dimension(SZIB_(HI),SZJ_(HI)), &
1706		optional, intent(inout) :: intx_dza !< The integral in x of the difference between
1707		!! the geopotential anomaly at the top and bottom of the layer divided
1708		!! by the x grid spacing [L2 T-2 ~> m2 s-2]
1709		real, dimension(SZI_(HI),SZJB_(HI)), &
1710		optional, intent(inout) :: inty_dza !< The integral in y of the difference between
1711		!! the geopotential anomaly at the top and bottom of the layer divided
1712		!! by the y grid spacing [L2 T-2 ~> m2 s-2]
1713		real, dimension(SZI_(HI),SZJ_(HI)), &
1714		optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
1715		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
1716		!! mass weighting to interpolate T/S in integrals
1717		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
1718		!! of T/S unless one side smaller than h_nv (i.e. vanished)
1719		real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]
1720
1721		! This subroutine calculates analytical and nearly-analytical integrals in
1722		! pressure across layers of geopotential anomalies, which are required for
1723		! calculating the finite-volume form pressure accelerations in a non-Boussinesq
1724		! model. There are essentially no free assumptions, apart from the use of
1725		! Boole's rule to do the horizontal integrals, and from a truncation in the
1726		! series for log(1-eps/1+eps) that assumes that \|eps\| < 0.34.
1727
1728	0	real :: T5((5HI%iscB+1):(5(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC]
1729	0	real :: S5((5HI%iscB+1):(5(HI%iecB+2))) ! Salinities along a line of subgrid locations [S ~> ppt]
1730	0	real :: p5((5HI%iscB+1):(5(HI%iecB+2))) ! Pressures along a line of subgrid locations [R L2 T-2 ~> Pa]
1731	0	real :: a5((5HI%iscB+1):(5(HI%iecB+2))) ! Specific volumes anomalies along a line of subgrid
1732		! locations [R-1 ~> m3 kg-1]
1733	0	real :: T15((15HI%iscB+1):(15(HI%iecB+1))) ! Temperatures at an array of subgrid locations [C ~> degC]
1734	0	real :: S15((15HI%iscB+1):(15(HI%iecB+1))) ! Salinities at an array of subgrid locations [S ~> ppt]
1735	0	real :: p15((15HI%iscB+1):(15(HI%iecB+1))) ! Pressures at an array of subgrid locations [R L2 T-2 ~> Pa]
1736	0	real :: a15((15HI%iscB+1):(15(HI%iecB+1))) ! Specific volumes at an array of subgrid locations [R ~> kg m-3]
1737		real :: wt_t(5), wt_b(5) ! Weights of top and bottom values at quadrature points [nondim]
1738		real :: T_top, T_bot ! Horizontally interpolated temperature at the cell top and bottom [C ~> degC]
1739		real :: S_top, S_bot ! Horizontally interpolated salinity at the cell top and bottom [S ~> ppt]
1740		real :: P_top, P_bot ! Horizontally interpolated pressure at the cell top and bottom [R L2 T-2 ~> Pa]
1741
1742		real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]
1743		real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]
1744	0	real :: dp_90(2:4,SZIB_(HI)) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa]
1745		real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
1746		real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
1747		real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
1748		real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim]
1749		real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]
1750		real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]
1751		real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]
1752		real :: intp(5) ! The integrals of specific volume with pressure at the
1753		! 5 sub-column locations [L2 T-2 ~> m2 s-2]
1754		real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
1755		logical :: do_massWeight ! Indicates whether to do mass weighting.
1756		logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
1757		logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting
1758		real :: massWeightNVonlyToggle ! A non-dimensional toggle factor for only using mass weighting
1759		! if at least one side vanished (0 or 1) [nondim]
1760		real :: p_nonvanished ! nonvanished pressure [R L2 T-2 ~> Pa]
1761		integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state
1762		integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state
1763		integer, dimension(2) :: EOSdom_h15 ! The 3x5-point h-point i-computational domain for the equation of state
1764		integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos
1765
1766	0	Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB
1767
1768	0	do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false.
1769	0	if (present(MassWghtInterp)) then
1770	0	do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
1771	0	top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
1772	0	massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set
1773	0	if (top_massWeight .and. .not.present(P_surf)) call MOM_error(FATAL, &
1774	0	"int_spec_vol_dp_generic_plm: P_surf must be present if near-surface mass weighting is in use.")
1775		endif
1776	0	massWeightNVonlyToggle = 1.
1777	0	if (present(MassWghtInterpVanOnly)) then
1778	0	if (MassWghtInterpVanOnly) massWeightNVonlyToggle = 0.
1779		endif
1780	0	p_nonvanished = 0.
1781	0	if (present(p_nv)) then
1782	0	p_nonvanished = p_nv
1783		endif
1784
1785	0	do n = 1, 5 ! Note that these are reversed from int_density_dz.
1786	0	wt_t(n) = 0.25 * real(n-1)
1787	0	wt_b(n) = 1.0 - wt_t(n)
1788		enddo
1789
1790		! Set the loop ranges for equation of state calculations at various points.
1791	0	EOSdom_h5(1) = 1 ; EOSdom_h5(2) = 5*(Ieq-Isq+2)
1792	0	EOSdom_q15(1) = 1 ; EOSdom_q15(2) = 15*(Ieq-Isq+1)
1793	0	EOSdom_h15(1) = 1 ; EOSdom_h15(2) = 15*(HI%iec-HI%isc+1)
1794
1795		! 1. Compute vertical integrals
1796	0	do j=Jsq,Jeq+1
1797	0	do i=Isq,Ieq+1
1798	0	do n=1,5 ! T, S and p are linearly interpolated in the vertical.
1799	0	p5(i5+n) = wt_t(n) p_t(i,j) + wt_b(n) * p_b(i,j)
1800	0	S5(i5+n) = wt_t(n) S_t(i,j) + wt_b(n) * S_b(i,j)
1801	0	T5(i5+n) = wt_t(n) T_t(i,j) + wt_b(n) * T_b(i,j)
1802		enddo
1803		enddo
1804	0	call calculate_spec_vol(T5, S5, p5, a5, EOS, EOSdom_h5, spv_ref=alpha_ref)
1805	0	do i=Isq,Ieq+1
1806		! Use Boole's rule to estimate the interface height anomaly change.
1807	0	dp = p_b(i,j) - p_t(i,j)
1808	0	alpha_anom = C1_90((7.0(a5(i5+1)+a5(i5+5)) + 32.0(a5(i5+2)+a5(i5+4))) + 12.0a5(i*5+3))
1809	0	dza(i,j) = dp*alpha_anom
1810		! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of
1811		! the interface height anomaly.
1812	0	if (present(intp_dza)) intp_dza(i,j) = 0.5dp2 &
1813	0	(alpha_anom - C1_90(16.0(a5(i5+4)-a5(i5+2)) + 7.0(a5(i5+5)-a5(i*5+1))) )
1814		enddo
1815		enddo
1816
1817		! 2. Compute horizontal integrals in the x direction
1818	0	if (present(intx_dza)) then ; do j=HI%jsc,HI%jec
1819	0	do I=Isq,Ieq
1820		! hWght is the distance measure by which the cell is violation of
1821		! hydrostatic consistency. For large hWght we bias the interpolation
1822		! of T,S along the top and bottom integrals, almost like thickness
1823		! weighting. Note: To work in terrain following coordinates we could
1824		! offset this distance by the layer thickness to replicate other models.
1825	0	hWght = 0.0
1826	0	if (do_massWeight .and. massWeight_bug) then
1827	0	hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j))
1828	0	elseif (do_massWeight) then
1829	0	hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j))
1830		endif
1831	0	if (top_massWeight) &
1832	0	hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j))
1833		! If both sides are nonvanished, then set it back to zero.
1834	0	if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i+1,j) - p_t(i+1,j)) > p_nonvanished)) then
1835	0	hWght = massWeightNVonlyToggle * hWght
1836		endif
1837	0	if (hWght > 0.) then
1838	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
1839	0	hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
1840	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
1841	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
1842	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
1843	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
1844		else
1845	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
1846		endif
1847
1848	0	do m=2,4
1849	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
1850	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
1851
1852		! T, S, and p are interpolated in the horizontal. The p interpolation
1853		! is linear, but for T and S it may be thickness weighted.
1854	0	P_top = (wt_Lp_t(i,j)) + (wt_Rp_t(i+1,j))
1855	0	P_bot = (wt_Lp_b(i,j)) + (wt_Rp_b(i+1,j))
1856	0	T_top = (wtT_LT_t(i,j)) + (wtT_RT_t(i+1,j))
1857	0	T_bot = (wtT_LT_b(i,j)) + (wtT_RT_b(i+1,j))
1858	0	S_top = (wtT_LS_t(i,j)) + (wtT_RS_t(i+1,j))
1859	0	S_bot = (wtT_LS_b(i,j)) + (wtT_RS_b(i+1,j))
1860	0	dp_90(m,I) = C1_90*(P_bot - P_top)
1861
1862		! Salinity, temperature and pressure with linear interpolation in the vertical.
1863	0	pos = i15+(m-2)5
1864	0	do n=1,5
1865	0	p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot
1866	0	S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot
1867	0	T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot
1868		enddo
1869		enddo
1870		enddo
1871
1872	0	call calculate_spec_vol(T15, S15, p15, a15, EOS, EOSdom_q15, spv_ref=alpha_ref)
1873
1874	0	do I=Isq,Ieq
1875	0	intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
1876	0	do m=2,4
1877		! Use Boole's rule to estimate the interface height anomaly change.
1878		! The integrals at the ends of the segment are already known.
1879	0	pos = I15+(m-2)5
1880		intp(m) = (dp_90(m,I)((7.0(a15(pos+1)+a15(pos+5)) + &
1881	0	32.0(a15(pos+2)+a15(pos+4))) + 12.0a15(pos+3) ))
1882		enddo
1883		! Use Boole's rule to integrate the interface height anomaly values in x.
1884		intx_dza(I,j) = C1_90((7.0(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
1885	0	12.0*intp(3))
1886		enddo
1887		enddo ; endif
1888
1889		! 3. Compute horizontal integrals in the y direction
1890	0	if (present(inty_dza)) then ; do J=Jsq,Jeq
1891	0	do i=HI%isc,HI%iec
1892		! hWght is the distance measure by which the cell is violation of
1893		! hydrostatic consistency. For large hWght we bias the interpolation
1894		! of T,S along the top and bottom integrals, like thickness weighting.
1895	0	hWght = 0.0
1896	0	if (do_massWeight .and. massWeight_bug) then
1897	0	hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j))
1898	0	elseif (do_massWeight) then
1899	0	hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1))
1900		endif
1901	0	if (top_massWeight) &
1902	0	hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j))
1903		! If both sides are nonvanished, then set it back to zero.
1904	0	if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i,j+1) - p_t(i,j+1)) > p_nonvanished)) then
1905	0	hWght = massWeightNVonlyToggle * hWght
1906		endif
1907	0	if (hWght > 0.) then
1908	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
1909	0	hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
1910	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
1911	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
1912	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
1913	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
1914		else
1915	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
1916		endif
1917
1918	0	do m=2,4
1919	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
1920	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
1921
1922		! T, S, and p are interpolated in the horizontal. The p interpolation
1923		! is linear, but for T and S it may be thickness weighted.
1924	0	P_top = (wt_Lp_t(i,j)) + (wt_Rp_t(i,j+1))
1925	0	P_bot = (wt_Lp_b(i,j)) + (wt_Rp_b(i,j+1))
1926	0	T_top = (wtT_LT_t(i,j)) + (wtT_RT_t(i,j+1))
1927	0	T_bot = (wtT_LT_b(i,j)) + (wtT_RT_b(i,j+1))
1928	0	S_top = (wtT_LS_t(i,j)) + (wtT_RS_t(i,j+1))
1929	0	S_bot = (wtT_LS_b(i,j)) + (wtT_RS_b(i,j+1))
1930	0	dp_90(m,i) = C1_90*(P_bot - P_top)
1931
1932		! Salinity, temperature and pressure with linear interpolation in the vertical.
1933	0	pos = i15+(m-2)5
1934	0	do n=1,5
1935	0	p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot
1936	0	S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot
1937	0	T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot
1938		enddo
1939		enddo
1940		enddo
1941
1942		call calculate_spec_vol(T15(15HI%isc+1:), S15(15HI%isc+1:), p15(15*HI%isc+1:), &
1943	0	a15(15*HI%isc+1:), EOS, EOSdom_h15, spv_ref=alpha_ref)
1944
1945	0	do i=HI%isc,HI%iec
1946	0	intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
1947	0	do m=2,4
1948		! Use Boole's rule to estimate the interface height anomaly change.
1949		! The integrals at the ends of the segment are already known.
1950	0	pos = i15+(m-2)5
1951		intp(m) = (dp_90(m,i) * ((7.0*(a15(pos+1)+a15(pos+5)) + &
1952	0	32.0(a15(pos+2)+a15(pos+4))) + 12.0a15(pos+3)))
1953		enddo
1954		! Use Boole's rule to integrate the interface height anomaly values in x.
1955		inty_dza(i,J) = C1_90((7.0(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
1956	0	12.0*intp(3))
1957		enddo
1958		enddo ; endif
1959
1960	0	end subroutine int_spec_vol_dp_generic_plm
1961
1962
1963		!> Diagnose the fractional mass weighting in a layer that might be used with a Boussinesq calculation.
1964	0	subroutine diagnose_mass_weight_Z(z_t, z_b, bathyT, SSH, dz_neglect, MassWghtInterp, HI, &
1965	0	MassWt_u, MassWt_v, MassWghtInterpVanOnly, h_nv)
1966		type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
1967		real, dimension(SZI_(HI),SZJ_(HI)), &
1968		intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]
1969		real, dimension(SZI_(HI),SZJ_(HI)), &
1970		intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]
1971		real, dimension(SZI_(HI),SZJ_(HI)), &
1972		intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]
1973		real, dimension(SZI_(HI),SZJ_(HI)), &
1974		intent(in) :: SSH !< The sea surface height [Z ~> m]
1975		real, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]
1976		integer, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
1977		!! mass weighting to interpolate T/S in integrals
1978		real, dimension(SZIB_(HI),SZJ_(HI)), &
1979		intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim]
1980		real, dimension(SZI_(HI),SZJB_(HI)), &
1981		intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim]
1982		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
1983		!! of T/S unless one side smaller than h_nv (i.e. vanished)
1984		real, optional, intent(in) :: h_nv !< Nonvanished height [Z ~> m]
1985
1986		! Local variables
1987		real :: hWght ! A pressure-thickness below topography [Z ~> m]
1988		real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m]
1989		real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2]
1990		logical :: do_massWeight ! Indicates whether to do mass weighting near bathymetry
1991		logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
1992		real :: massWeightNVonlyToggle ! A non-dimensional toggle factor for only using mass weighting
1993		real :: h_nonvanished ! nonvanished height [Z ~> m]
1994		integer :: Isq, Ieq, Jsq, Jeq, i, j
1995
1996	0	Isq = HI%IscB ; Ieq = HI%IecB
1997	0	Jsq = HI%JscB ; Jeq = HI%JecB
1998
1999	0	do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
2000	0	top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
2001	0	massWeightNVonlyToggle = 1.
2002	0	if (present(MassWghtInterpVanOnly)) then
2003	0	if (MassWghtInterpVanOnly) massWeightNVonlyToggle = 0.
2004		endif
2005	0	h_nonvanished = 0.
2006	0	if (present(h_nv)) then
2007	0	h_nonvanished = h_nv
2008		endif
2009
2010		! Calculate MassWt_u
2011	0	do j=HI%jsc,HI%jec ; do I=Isq,Ieq
2012		! hWght is the distance measure by which the cell is violation of
2013		! hydrostatic consistency. For large hWght we bias the interpolation
2014		! of T,S along the top and bottom integrals, like thickness weighting.
2015	0	hWght = 0.0
2016	0	if (do_massWeight) &
2017	0	hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j))
2018	0	if (top_massWeight) &
2019	0	hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j))
2020		! If both sides are nonvanished, then set it back to zero.
2021	0	if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i+1,j) - z_b(i+1,j)) > h_nonvanished)) then
2022	0	hWght = massWeightNVonlyToggle * hWght
2023		endif
2024	0	if (hWght > 0.) then
2025	0	hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
2026	0	hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
2027	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
2028	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
2029	0	MassWt_u(I,j) = (hWghthR + hWghthL) * iDenom
2030		else
2031	0	MassWt_u(I,j) = 0.0
2032		endif
2033		enddo ; enddo
2034
2035		! Calculate MassWt_v
2036	0	do J=Jsq,Jeq ; do i=HI%isc,HI%iec
2037		! hWght is the distance measure by which the cell is violation of
2038		! hydrostatic consistency. For large hWght we bias the interpolation
2039		! of T,S along the top and bottom integrals, like thickness weighting.
2040	0	hWght = 0.0
2041	0	if (do_massWeight) &
2042	0	hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j))
2043	0	if (top_massWeight) &
2044	0	hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1))
2045		! If both sides are nonvanished, then set it back to zero.
2046	0	if (((z_t(i,j) - z_b(i,j)) > h_nonvanished) .and. ((z_t(i,j+1) - z_b(i,j+1)) > h_nonvanished)) then
2047	0	hWght = massWeightNVonlyToggle * hWght
2048		endif
2049	0	if (hWght > 0.) then
2050	0	hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
2051	0	hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
2052	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
2053	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
2054	0	MassWt_v(i,J) = (hWghthR + hWghthL) * iDenom
2055		else
2056	0	MassWt_v(i,J) = 0.0
2057		endif
2058		enddo ; enddo
2059
2060	0	end subroutine diagnose_mass_weight_Z
2061
2062
2063		!> Diagnose the fractional mass weighting in a layer that might be used with a non-Boussinesq calculation.
2064	0	subroutine diagnose_mass_weight_p(p_t, p_b, bathyP, P_surf, dP_neglect, MassWghtInterp, HI, &
2065	0	MassWt_u, MassWt_v, MassWghtInterpVanOnly, p_nv)
2066		type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
2067		real, dimension(SZI_(HI),SZJ_(HI)), &
2068		intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa]
2069		real, dimension(SZI_(HI),SZJ_(HI)), &
2070		intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa]
2071		real, intent(in) :: dP_neglect !<!< A miniscule pressure change with
2072		!! the same units as p_t [R L2 T-2 ~> Pa]
2073		real, dimension(SZI_(HI),SZJ_(HI)), &
2074		intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
2075		real, dimension(SZI_(HI),SZJ_(HI)), &
2076		intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
2077		integer, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
2078		!! mass weighting to interpolate T/S in integrals
2079		real, dimension(SZIB_(HI),SZJ_(HI)), &
2080		intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim]
2081		real, dimension(SZI_(HI),SZJB_(HI)), &
2082		intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim]
2083		logical, optional, intent(in) :: MassWghtInterpVanOnly !< If true, does not do mass weighting
2084		!! of T/S unless one side smaller than h_nv (i.e. vanished)
2085		real, optional, intent(in) :: p_nv !< Nonvanished pressure [R L2 T-2 ~> Pa]
2086
2087		! Local variables
2088		real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
2089		real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
2090		real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
2091		logical :: do_massWeight ! Indicates whether to do mass weighting.
2092		logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
2093		logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting
2094		real :: massWeightNVonlyToggle ! A non-dimensional toggle factor for only using mass weighting
2095		! if at least one side vanished (0 or 1) [nondim]
2096		real :: p_nonvanished ! nonvanished pressure [R L2 T-2 ~> Pa]
2097
2098		integer :: Isq, Ieq, Jsq, Jeq, i, j
2099
2100	0	Isq = HI%IscB ; Ieq = HI%IecB
2101	0	Jsq = HI%JscB ; Jeq = HI%JecB
2102
2103	0	do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
2104	0	top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
2105	0	massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set
2106	0	massWeightNVonlyToggle = 1.
2107	0	if (present(MassWghtInterpVanOnly)) then
2108	0	if (MassWghtInterpVanOnly) massWeightNVonlyToggle = 0.
2109		endif
2110	0	p_nonvanished = 0.
2111	0	if (present(p_nv)) then
2112	0	p_nonvanished = p_nv
2113		endif
2114
2115		! Calculate MassWt_u
2116	0	do j=HI%jsc,HI%jec ; do I=Isq,Ieq
2117		! hWght is the distance measure by which the cell is violation of
2118		! hydrostatic consistency. For large hWght we bias the interpolation
2119		! of T,S along the top and bottom integrals, like thickness weighting.
2120	0	hWght = 0.0
2121	0	if (do_massWeight .and. massWeight_bug) then
2122	0	hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j))
2123	0	elseif (do_massWeight) then
2124	0	hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j))
2125		endif
2126	0	if (top_massWeight) &
2127	0	hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j))
2128		! If both sides are nonvanished, then set it back to zero.
2129	0	if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i+1,j) - p_t(i+1,j)) > p_nonvanished)) then
2130	0	hWght = massWeightNVonlyToggle * hWght
2131		endif
2132	0	if (hWght > 0.) then
2133	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
2134	0	hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
2135	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
2136	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
2137	0	MassWt_u(I,j) = (hWghthR + hWghthL) * iDenom
2138		else
2139	0	MassWt_u(I,j) = 0.0
2140		endif
2141		enddo ; enddo
2142
2143		! Calculate MassWt_v
2144	0	do J=Jsq,Jeq ; do i=HI%isc,HI%iec
2145		! hWght is the distance measure by which the cell is violation of
2146		! hydrostatic consistency. For large hWght we bias the interpolation
2147		! of T,S along the top and bottom integrals, like thickness weighting.
2148	0	hWght = 0.0
2149	0	if (do_massWeight .and. massWeight_bug) then
2150	0	hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j))
2151	0	elseif (do_massWeight) then
2152	0	hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1))
2153		endif
2154	0	if (top_massWeight) &
2155	0	hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j))
2156		! If both sides are nonvanished, then set it back to zero.
2157	0	if (((p_b(i,j) - p_t(i,j)) > p_nonvanished) .and. ((p_b(i,j+1) - p_t(i,j+1)) > p_nonvanished)) then
2158	0	hWght = massWeightNVonlyToggle * hWght
2159		endif
2160	0	if (hWght > 0.) then
2161	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
2162	0	hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
2163	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
2164	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
2165	0	MassWt_v(i,J) = (hWghthR + hWghthL) * iDenom
2166		else
2167	0	MassWt_v(i,J) = 0.0
2168		endif
2169		enddo ; enddo
2170
2171	0	end subroutine diagnose_mass_weight_p
2172
2173		!> Find the depth at which the reconstructed pressure matches P_tgt
2174	0	subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, &
2175		rho_ref, G_e, EOS, US, P_b, z_out, z_tol, frac_dp_bugfix)
2176		real, intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC]
2177		real, intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC]
2178		real, intent(in) :: S_t !< Salinity at the cell top [S ~> ppt]
2179		real, intent(in) :: S_b !< Salinity at the cell bottom [S ~> ppt]
2180		real, intent(in) :: z_t !< Absolute height of top of cell [Z ~> m] (Boussinesq ????)
2181		real, intent(in) :: z_b !< Absolute height of bottom of cell [Z ~> m]
2182		real, intent(in) :: P_t !< Anomalous pressure of top of cell, relative
2183		!! to grho_refz_t [R L2 T-2 ~> Pa]
2184		real, intent(in) :: P_tgt !< Target pressure at height z_out, relative
2185		!! to grho_refz_out [R L2 T-2 ~> Pa]
2186		real, intent(in) :: rho_ref !< Reference density with which calculation
2187		!! are anomalous to [R ~> kg m-3]
2188		real, intent(in) :: G_e !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
2189		type(EOS_type), intent(in) :: EOS !< Equation of state structure
2190		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2191		real, intent(out) :: P_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa]
2192		real, intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m]
2193		real, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m]
2194		logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
2195
2196		! Local variables
2197		real :: dp ! Pressure thickness of the layer [R L2 T-2 ~> Pa]
2198		real :: F_guess, F_l, F_r ! Fractional positions [nondim]
2199		real :: GxRho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]
2200		real :: Pa, Pa_left, Pa_right, Pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz [R L2 T-2 ~> Pa]
2201		integer :: m ! A counter for how many iterations have been done in the while loop
2202		character(len=240) :: msg
2203
2204	0	GxRho = G_e * rho_ref
2205
2206		! Anomalous pressure difference across whole cell
2207	0	dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS, frac_dp_bugfix)
2208
2209	0	P_b = P_t + dp ! Anomalous pressure at bottom of cell
2210
2211	0	if (P_tgt <= P_t ) then
2212	0	z_out = z_t
2213	0	return
2214		endif
2215
2216	0	if (P_tgt >= P_b) then
2217	0	z_out = z_b
2218	0	return
2219		endif
2220
2221	0	F_l = 0.
2222	0	Pa_left = P_t - P_tgt ! Pa_left < 0
2223	0	F_r = 1.
2224	0	Pa_right = P_b - P_tgt ! Pa_right > 0
2225	0	Pa_tol = GxRho * z_tol
2226
2227	0	F_guess = F_l - Pa_left / (Pa_right - Pa_left) * (F_r - F_l)
2228	0	Pa = Pa_right - Pa_left ! To get into iterative loop
2229	0	m = 0 ! Reset the counter for the loop to be zero
2230	0	do while ( abs(Pa) > Pa_tol )
2231
2232	0	m = m + 1
2233	0	if (m > 30) then ! Call an error, because convergence to the tolerance has not been achieved
2234	0	write(msg,*) Pa_left,Pa,Pa_right,P_t-P_tgt,P_b-P_tgt
2235	0	call MOM_error(FATAL, 'find_depth_of_pressure_in_cell completes too many iterations: '//msg)
2236		endif
2237	0	z_out = z_t + ( z_b - z_t ) * F_guess
2238	0	Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS, frac_dp_bugfix) - ( P_tgt - P_t )
2239
2240	0	if (Pa<Pa_left) then
2241	0	write(msg,*) Pa_left,Pa,Pa_right,P_t-P_tgt,P_b-P_tgt
2242	0	call MOM_error(FATAL, 'find_depth_of_pressure_in_cell out of bounds negative: /n'//msg)
2243	0	elseif (Pa<0.) then
2244	0	Pa_left = Pa
2245	0	F_l = F_guess
2246	0	elseif (Pa>Pa_right) then
2247	0	write(msg,*) Pa_left,Pa,Pa_right,P_t-P_tgt,P_b-P_tgt
2248	0	call MOM_error(FATAL, 'find_depth_of_pressure_in_cell out of bounds positive: /n'//msg)
2249	0	elseif (Pa>0.) then
2250	0	Pa_right = Pa
2251	0	F_r = F_guess
2252		else ! Pa == 0
2253	0	return
2254		endif
2255	0	F_guess = F_l - Pa_left / (Pa_right - Pa_left) * (F_r - F_l)
2256
2257		enddo
2258
2259		end subroutine find_depth_of_pressure_in_cell
2260
2261		!> Calculate the average in situ specific volume across layers
2262	0	subroutine avg_specific_vol(T, S, p_t, dp, HI, EOS, SpV_avg, halo_size)
2263		type(hor_index_type), intent(in) :: HI !< The horizontal index structure
2264		real, dimension(SZI_(HI),SZJ_(HI)), &
2265		intent(in) :: T !< Potential temperature of the layer [C ~> degC]
2266		real, dimension(SZI_(HI),SZJ_(HI)), &
2267		intent(in) :: S !< Salinity of the layer [S ~> ppt]
2268		real, dimension(SZI_(HI),SZJ_(HI)), &
2269		intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
2270		real, dimension(SZI_(HI),SZJ_(HI)), &
2271		intent(in) :: dp !< Pressure change in the layer [R L2 T-2 ~> Pa]
2272		type(EOS_type), intent(in) :: EOS !< Equation of state structure
2273		real, dimension(SZI_(HI),SZJ_(HI)), &
2274		intent(inout) :: SpV_avg !< The vertical average specific volume
2275		!! in the layer [R-1 ~> m3 kg-1]
2276		integer, optional, intent(in) :: halo_size !< The number of halo points in which to work.
2277
2278		! Local variables
2279		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2280		integer :: jsh, jeh, j, halo
2281
2282	0	halo = 0 ; if (present(halo_size)) halo = MAX(halo_size,0)
2283	0	jsh = HI%jsc-halo ; jeh = HI%jec+halo
2284
2285	0	EOSdom(:) = EOS_domain(HI, halo_size)
2286	0	do j=jsh,jeh
2287	0	call average_specific_vol(T(:,j), S(:,j), p_t(:,j), dp(:,j), SpV_avg(:,j), EOS, EOSdom)
2288		enddo
2289
2290	0	end subroutine avg_specific_vol
2291
2292		!> Returns change in anomalous pressure change from top to non-dimensional
2293		!! position pos between z_t and z_b [R L2 T-2 ~> Pa]
2294	0	real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS, frac_dp_bugfix)
2295		real, intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC]
2296		real, intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC]
2297		real, intent(in) :: S_t !< Salinity at the cell top [S ~> ppt]
2298		real, intent(in) :: S_b !< Salinity at the cell bottom [S ~> ppt]
2299		real, intent(in) :: z_t !< The geometric height at the top of the layer [Z ~> m]
2300		real, intent(in) :: z_b !< The geometric height at the bottom of the layer [Z ~> m]
2301		real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted out to
2302		!! reduce the magnitude of each of the integrals.
2303		real, intent(in) :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
2304		real, intent(in) :: pos !< The fractional vertical position, 0 to 1 [nondim]
2305		type(EOS_type), intent(in) :: EOS !< Equation of state structure
2306		logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
2307
2308		! Local variables
2309		real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
2310		real :: dz ! Distance from the layer top [Z ~> m]
2311		real :: top_weight, bottom_weight ! Fractional weights at quadrature points [nondim]
2312		real :: rho_ave ! Average density [R ~> kg m-3]
2313		real, dimension(5) :: T5 ! Temperatures at quadrature points [C ~> degC]
2314		real, dimension(5) :: S5 ! Salinities at quadrature points [S ~> ppt]
2315		real, dimension(5) :: p5 ! Pressures at quadrature points [R L2 T-2 ~> Pa]
2316		real, dimension(5) :: rho5 ! Densities at quadrature points [R ~> kg m-3]
2317		integer :: n
2318
2319	0	do n=1,5
2320		! Evaluate density at five quadrature points
2321	0	bottom_weight = 0.25real(n-1) pos
2322	0	top_weight = 1.0 - bottom_weight
2323		! Salinity and temperature points are linearly interpolated
2324	0	S5(n) = top_weight * S_t + bottom_weight * S_b
2325	0	T5(n) = top_weight * T_t + bottom_weight * T_b
2326	0	if (frac_dp_bugfix) then
2327	0	p5(n) = (-1) * ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
2328		else
2329	0	p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
2330		endif !bugfix
2331		enddo
2332	0	call calculate_density(T5, S5, p5, rho5, EOS)
2333	0	rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref
2334
2335		! Use Boole's rule to estimate the average density
2336	0	rho_ave = C1_90(7.0(rho5(1)+rho5(5)) + 32.0(rho5(2)+rho5(4)) + 12.0rho5(3))
2337
2338	0	dz = ( z_t - z_b ) * pos
2339	0	frac_dp_at_pos = G_e * dz * rho_ave
2340	0	end function frac_dp_at_pos
2341
2342		end module MOM_density_integrals
2343
2344		!> \namespace mom_density_integrals
2345		!!