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