src/equation_of_state/MOM_EOS_Wright.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		!> The equation of state using a poor implementation (missing parenthesis and bugs) of the
6		!! reduced range Wright 1997 expressions
7		module MOM_EOS_Wright
8
9		use MOM_EOS_base_type, only : EOS_base
10		use MOM_hor_index, only : hor_index_type
11
12		implicit none ; private
13
14		public buggy_Wright_EOS
15		public int_density_dz_wright, int_spec_vol_dp_wright
16		public avg_spec_vol_buggy_Wright
17		public set_params_buggy_Wright
18
19		!>@{ Parameters in the Wright equation of state using the reduced range formula, which is a fit to the UNESCO
20		! equation of state for the restricted range: -2 < theta < 30 [degC], 28 < S < 38 [PSU], 0 < p < 5e7 [Pa].
21
22		! Note that a0/a1 ~= 2028 [degC] ; a0/a2 ~= -6343 [PSU]
23		! b0/b1 ~= 165 [degC] ; b0/b4 ~= 974 [PSU]
24		! c0/c1 ~= 216 [degC] ; c0/c4 ~= -740 [PSU]
25		real, parameter :: a0 = 7.057924e-4 ! A parameter in the Wright alpha_0 fit [m3 kg-1]
26		real, parameter :: a1 = 3.480336e-7 ! A parameter in the Wright alpha_0 fit [m3 kg-1 degC-1]
27		real, parameter :: a2 = -1.112733e-7 ! A parameter in the Wright alpha_0 fit [m3 kg-1 PSU-1]
28		real, parameter :: b0 = 5.790749e8 ! A parameter in the Wright p_0 fit [Pa]
29		real, parameter :: b1 = 3.516535e6 ! A parameter in the Wright p_0 fit [Pa degC-1]
30		real, parameter :: b2 = -4.002714e4 ! A parameter in the Wright p_0 fit [Pa degC-2]
31		real, parameter :: b3 = 2.084372e2 ! A parameter in the Wright p_0 fit [Pa degC-3]
32		real, parameter :: b4 = 5.944068e5 ! A parameter in the Wright p_0 fit [Pa PSU-1]
33		real, parameter :: b5 = -9.643486e3 ! A parameter in the Wright p_0 fit [Pa degC-1 PSU-1]
34		real, parameter :: c0 = 1.704853e5 ! A parameter in the Wright lambda fit [m2 s-2]
35		real, parameter :: c1 = 7.904722e2 ! A parameter in the Wright lambda fit [m2 s-2 degC-1]
36		real, parameter :: c2 = -7.984422 ! A parameter in the Wright lambda fit [m2 s-2 degC-2]
37		real, parameter :: c3 = 5.140652e-2 ! A parameter in the Wright lambda fit [m2 s-2 degC-3]
38		real, parameter :: c4 = -2.302158e2 ! A parameter in the Wright lambda fit [m2 s-2 PSU-1]
39		real, parameter :: c5 = -3.079464 ! A parameter in the Wright lambda fit [m2 s-2 degC-1 PSU-1]
40		!>@}
41
42		!> The EOS_base implementation of the Wright 1997 equation of state with some bugs
43		type, extends (EOS_base) :: buggy_Wright_EOS
44
45		real :: three = 3.0 !< A constant that can be adjusted to recreate some bugs [nondim]
46
47		contains
48		!> Implementation of the in-situ density as an elemental function [kg m-3]
49		procedure :: density_elem => density_elem_buggy_Wright
50		!> Implementation of the in-situ density anomaly as an elemental function [kg m-3]
51		procedure :: density_anomaly_elem => density_anomaly_elem_buggy_Wright
52		!> Implementation of the in-situ specific volume as an elemental function [m3 kg-1]
53		procedure :: spec_vol_elem => spec_vol_elem_buggy_Wright
54		!> Implementation of the in-situ specific volume anomaly as an elemental function [m3 kg-1]
55		procedure :: spec_vol_anomaly_elem => spec_vol_anomaly_elem_buggy_Wright
56		!> Implementation of the calculation of derivatives of density
57		procedure :: calculate_density_derivs_elem => calculate_density_derivs_elem_buggy_Wright
58		!> Implementation of the calculation of second derivatives of density
59		procedure :: calculate_density_second_derivs_elem => calculate_density_second_derivs_elem_buggy_Wright
60		!> Implementation of the calculation of derivatives of specific volume
61		procedure :: calculate_specvol_derivs_elem => calculate_specvol_derivs_elem_buggy_Wright
62		!> Implementation of the calculation of compressibility
63		procedure :: calculate_compress_elem => calculate_compress_elem_buggy_Wright
64		!> Implementation of the range query function
65		procedure :: EOS_fit_range => EOS_fit_range_buggy_Wright
66
67		!> Instance specific function to set internal parameters
68		procedure :: set_params_buggy_Wright => set_params_buggy_Wright
69
70		!> Local implementation of generic calculate_density_array for efficiency
71		procedure :: calculate_density_array => calculate_density_array_buggy_Wright
72		!> Local implementation of generic calculate_density_array_2d for efficiency
73		procedure :: calculate_density_array_2d => calculate_density_array_2d_buggy_Wright
74		!> Local implementation of generic calculate_density_array_3d for efficiency
75		procedure :: calculate_density_array_3d => calculate_density_array_3d_buggy_Wright
76		!> Local implementation of generic calculate_spec_vol_array for efficiency
77		procedure :: calculate_spec_vol_array => calculate_spec_vol_array_buggy_Wright
78		!> Local implementation of generic calculate_density_derivs_2d for efficiency
79		procedure :: calculate_density_derivs_2d => calculate_density_derivs_2d_buggy_Wright
80		!> Local implementation of generic calculate_density_derivs_3d for efficiency
81		procedure :: calculate_density_derivs_3d => calculate_density_derivs_3d_buggy_Wright
82
83		end type buggy_Wright_EOS
84
85		contains
86
87		!> In situ density of sea water using a buggy implementation of Wright, 1997 [kg m-3]
88		!!
89		!! This is an elemental function that can be applied to any combination of scalar and array inputs.
90	7065260	real elemental function density_elem_buggy_Wright_loc(T, S, pressure)
91		real, intent(in) :: T !< potential temperature relative to the surface [degC].
92		real, intent(in) :: S !< salinity [PSU].
93		real, intent(in) :: pressure !< pressure [Pa].
94
95		! Local variables
96		real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
97		real :: p0 ! The pressure offset in the Wright EOS [Pa]
98		real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
99
100	7065260	al0 = (a0 + a1T) +a2S
101	7065260	p0 = (b0 + b4S) + T (b1 + T(b2 + b3T) + b5*S)
102	7065260	lambda = (c0 +c4S) + T (c1 + T(c2 + c3T) + c5*S)
103	7065260	density_elem_buggy_Wright_loc = (pressure + p0) / (lambda + al0*(pressure + p0))
104
105	7065260	end function density_elem_buggy_Wright_loc
106
107		!> Wrapper for density_elem_buggy_Wright_loc created to preserve API while calling
108		!! density_elem_buggy_Wright without "this" variable that causes runtime errors on
109		!! gpu runs with nvfortran.
110	2	real elemental function density_elem_buggy_Wright(this, T, S, pressure)
111		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
112		real, intent(in) :: T !< potential temperature relative to the surface [degC].
113		real, intent(in) :: S !< salinity [PSU].
114		real, intent(in) :: pressure !< pressure [Pa].
115
116	2	density_elem_buggy_Wright = density_elem_buggy_Wright_loc(T, S, pressure)
117	2	end function density_elem_buggy_Wright
118
119		!> In situ density anomaly of sea water using a buggy implementation of Wright, 1997 [kg m-3]
120		!!
121		!! This is an elemental function that can be applied to any combination of scalar and array inputs.
122	0	real elemental function density_anomaly_elem_buggy_Wright(this, T, S, pressure, rho_ref)
123		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
124		real, intent(in) :: T !< potential temperature relative to the surface [degC].
125		real, intent(in) :: S !< salinity [PSU].
126		real, intent(in) :: pressure !< pressure [Pa].
127		real, intent(in) :: rho_ref !< A reference density [kg m-3].
128
129		! Local variables
130		real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
131		real :: p0 ! The pressure offset in the Wright EOS [Pa]
132		real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
133		real :: al_TS ! The contributions of temperature and salinity to al0 [m3 kg-1]
134		real :: p_TSp ! A combination of the pressure and the temperature and salinity contributions to p0 [Pa]
135		real :: lam_TS ! The contributions of temperature and salinity to lambda [m2 s-2]
136		real :: pa_000 ! A corrected offset to the pressure, including contributions from rho_ref [Pa]
137
138	0	pa_000 = (b0(1.0 - a0rho_ref) - rho_ref*c0)
139	0	al_TS = a1T +a2S
140	0	al0 = a0 + al_TS
141	0	p_TSp = pressure + (b4S + T (b1 + (T(b2 + b3T) + b5*S)))
142	0	lam_TS = c4S + T (c1 + (T(c2 + c3T) + c5*S))
143
144		! The following two expressions are mathematically equivalent.
145		! wright_density = (b0 + p0_TSp) / ((c0 + lam_TS) + al0*(b0 + p0_TSp)) - rho_ref
146		density_anomaly_elem_buggy_Wright = &
147		(pa_000 + (p_TSp - rho_ref(p_TSpal0 + (b0*al_TS + lam_TS)))) / &
148	0	( (c0 + lam_TS) + al0*(b0 + p_TSp) )
149
150	0	end function density_anomaly_elem_buggy_Wright
151
152		!> In situ specific volume of sea water using a buggy implementation of Wright, 1997 [kg m-3]
153		!!
154		!! This is an elemental function that can be applied to any combination of scalar and array inputs.
155	0	real elemental function spec_vol_elem_buggy_Wright(this, T, S, pressure)
156		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
157		real, intent(in) :: T !< potential temperature relative to the surface [degC].
158		real, intent(in) :: S !< salinity [PSU].
159		real, intent(in) :: pressure !< pressure [Pa].
160
161		! Local variables
162		real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
163		real :: p0 ! The pressure offset in the Wright EOS [Pa]
164		real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
165
166	0	al0 = (a0 + a1T) +a2S
167	0	p0 = (b0 + b4S) + T (b1 + T((b2 + b3T)) + b5*S)
168	0	lambda = (c0 +c4S) + T (c1 + T((c2 + c3T)) + c5*S)
169
170	0	spec_vol_elem_buggy_Wright = (lambda + al0*(pressure + p0)) / (pressure + p0)
171
172	0	end function spec_vol_elem_buggy_Wright
173
174		!> In situ specific volume anomaly of sea water using a buggy implementation of Wright, 1997 [kg m-3]
175		!!
176		!! This is an elemental function that can be applied to any combination of scalar and array inputs.
177	0	real elemental function spec_vol_anomaly_elem_buggy_Wright(this, T, S, pressure, spv_ref)
178		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
179		real, intent(in) :: T !< potential temperature relative to the surface [degC].
180		real, intent(in) :: S !< salinity [PSU].
181		real, intent(in) :: pressure !< pressure [Pa].
182		real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
183
184		! Local variables
185		real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
186		real :: p0 ! The pressure offset in the Wright EOS [Pa]
187		real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
188
189	0	al0 = (a0 + a1T) +a2S
190	0	p0 = (b0 + b4S) + T (b1 + T((b2 + b3T)) + b5*S)
191	0	lambda = (c0 +c4S) + T (c1 + T((c2 + c3T)) + c5*S)
192
193	0	spec_vol_anomaly_elem_buggy_Wright = (lambda + (al0 - spv_ref)*(pressure + p0)) / (pressure + p0)
194
195	0	end function spec_vol_anomaly_elem_buggy_Wright
196
197		!> Calculate the partial derivatives of density with potential temperature and salinity
198		!! using the buggy implementation of the equation of state, as fit by Wright, 1997
199	10966107	elemental subroutine calculate_density_derivs_elem_buggy_Wright_loc( T, S, pressure, drho_dT, drho_dS)
200		real, intent(in) :: T !< Potential temperature relative to the surface [degC]
201		real, intent(in) :: S !< Salinity [PSU]
202		real, intent(in) :: pressure !< Pressure [Pa]
203		real, intent(out) :: drho_dT !< The partial derivative of density with potential
204		!! temperature [kg m-3 degC-1]
205		real, intent(out) :: drho_dS !< The partial derivative of density with salinity,
206		!! in [kg m-3 PSU-1]
207		! Local variables
208		real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
209		real :: p0 ! The pressure offset in the Wright EOS [Pa]
210		real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
211		real :: I_denom2 ! The inverse of the square of the denominator of density in the Wright EOS [s4 m-4]
212
213	10966107	al0 = (a0 + a1T) + a2S
214	10966107	p0 = (b0 + b4S) + T (b1 + T((b2 + b3T)) + b5*S)
215	10966107	lambda = (c0 +c4S) + T (c1 + T((c2 + c3T)) + c5*S)
216
217	10966107	I_denom2 = 1.0 / (lambda + al0*(pressure + p0))
218	10966107	I_denom2 = I_denom2 *I_denom2
219		drho_dT = I_denom2 * &
220		(lambda* (b1 + T(2.0b2 + 3.0b3T) + b5*S) - &
221		(pressure+p0) * ( (pressure+p0)*a1 + &
222	10966107	(c1 + T(c22.0 + c33.0T) + c5*S) ))
223		drho_dS = I_denom2 * (lambda* (b4 + b5*T) - &
224	10966107	(pressure+p0) * ( (pressure+p0)a2 + (c4 + c5T) ))
225
226	10966107	end subroutine calculate_density_derivs_elem_buggy_Wright_loc
227
228		!> Wrapper for density_elem_buggy_Wright_loc created to preserve API while calling
229		!! density_elem_buggy_Wright without "this" variable that causes runtime errors on
230		!! gpu runs with nvfortran.
231	5431851	elemental subroutine calculate_density_derivs_elem_buggy_Wright(this, T, S, pressure, drho_dT, drho_dS)
232		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
233		real, intent(in) :: T !< Potential temperature relative to the surface [degC]
234		real, intent(in) :: S !< Salinity [PSU]
235		real, intent(in) :: pressure !< Pressure [Pa]
236		real, intent(out) :: drho_dT !< The partial derivative of density with potential
237		!! temperature [kg m-3 degC-1]
238		real, intent(out) :: drho_dS !< The partial derivative of density with salinity,
239		!! in [kg m-3 PSU-1]
240
241	5431851	call calculate_density_derivs_elem_buggy_Wright_loc(T, S, pressure, drho_dT, drho_dS)
242
243	5431851	end subroutine calculate_density_derivs_elem_buggy_Wright
244
245		!> Second derivatives of density with respect to temperature, salinity, and pressure,
246		!! using the poor implementation of the equation of state, as fit by Wright, 1997
247	0	elemental subroutine calculate_density_second_derivs_elem_buggy_Wright(this, T, S, pressure, &
248		drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
249		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
250		real, intent(in) :: T !< Potential temperature referenced to 0 dbar [degC]
251		real, intent(in) :: S !< Salinity [PSU]
252		real, intent(in) :: pressure !< Pressure [Pa]
253		real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
254		!! to S [kg m-3 PSU-2]
255		real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
256		!! to T [kg m-3 PSU-1 degC-1]
257		real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
258		!! to T [kg m-3 degC-2]
259		real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
260		!! to pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1]
261		real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
262		!! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
263		! Local variables
264		real :: z0, z1 ! Local work variables [Pa]
265		real :: z2, z4 ! Local work variables [m2 s-2]
266		real :: z3, z5 ! Local work variables [Pa degC-1]
267		real :: z6, z8 ! Local work variables [m2 s-2 degC-1]
268		real :: z7 ! A local work variable [m2 s-2 PSU-1]
269		real :: z9 ! A local work variable [m3 kg-1]
270		real :: z10 ! A local work variable [Pa PSU-1]
271		real :: z11 ! A local work variable [Pa m2 s-2 PSU-1] = [kg m s-4 PSU-1]
272		real :: z2_2 ! A local work variable [m4 s-4]
273		real :: z2_3 ! A local work variable [m6 s-6]
274		real :: six ! A constant that can be adjusted from 6. to 4. to recreate a bug [nondim]
275
276		! Based on the above expression with common terms factored, there probably exists a more numerically stable
277		! and/or efficient expression
278
279	0	six = 2.0*this%three ! When recreating a bug from the original version of this routine, six = 4.
280
281	0	z0 = T(b1 + b5S + T(b2 + b3T))
282	0	z1 = (b0 + pressure + b4*S + z0)
283	0	z3 = (b1 + b5S + T(2.b2 + this%threeb3*T)) ! When recreating a bug here this%three = 2.
284	0	z4 = (c0 + c4S + T(c1 + c5S + T(c2 + c3*T)))
285	0	z5 = (b1 + b5S + T(b2 + b3T) + T(b2 + 2.b3T))
286	0	z6 = c1 + c5S + T(c2 + c3T) + T(c2 + 2.c3T)
287	0	z7 = (c4 + c5T + a2z1)
288	0	z8 = (c1 + c5S + T(2.c2 + 3.c3T) + a1z1)
289	0	z9 = (a0 + a2S + a1T)
290	0	z10 = (b4 + b5*T)
291	0	z11 = (z10z4 - z1z7)
292	0	z2 = (c0 + c4S + T(c1 + c5S + T(c2 + c3T)) + z9z1)
293	0	z2_2 = z2*z2
294	0	z2_3 = z2_2*z2
295
296	0	drho_ds_ds = (z10(c4 + c5T) - a2z10z1 - z10z7)/z2_2 - (2.(c4 + c5T + z9z10 + a2z1)z11)/z2_3
297	0	drho_ds_dt = (z10z6 - z1(c5 + a2z5) + b5z4 - z5z7)/z2_2 - (2.(z6 + z9z5 + a1z1)*z11)/z2_3
298		drho_dt_dt = (z3z6 - z1(2.c2 + 6.c3T + a1z5) + (2.b2 + sixb3T)z4 - z5*z8)/z2_2 - &
299	0	(2.(z6 + z9z5 + a1z1)(z3z4 - z1z8))/z2_3
300	0	drho_ds_dp = (-c4 - c5T - 2.a2z1)/z2_2 - (2.z9*z11)/z2_3
301	0	drho_dt_dp = (-c1 - c5S - T(2.c2 + 3.c3T) - 2.a1z1)/z2_2 - (2.z9(z3z4 - z1*z8))/z2_3
302
303	0	end subroutine calculate_density_second_derivs_elem_buggy_Wright
304
305		!> Calculate the partial derivatives of specific volume with temperature and salinity
306		!! using the poor implementation of the equation of state, as fit by Wright, 1997
307	0	elemental subroutine calculate_specvol_derivs_elem_buggy_Wright(this, T, S, pressure, dSV_dT, dSV_dS)
308		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
309		real, intent(in) :: T !< Potential temperature [degC]
310		real, intent(in) :: S !< Salinity [PSU]
311		real, intent(in) :: pressure !< Pressure [Pa]
312		real, intent(inout) :: dSV_dT !< The partial derivative of specific volume with
313		!! potential temperature [m3 kg-1 degC-1]
314		real, intent(inout) :: dSV_dS !< The partial derivative of specific volume with
315		!! salinity [m3 kg-1 PSU-1]
316		! Local variables
317		real :: p0 ! The pressure offset in the Wright EOS [Pa]
318		real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
319		real :: I_denom ! The inverse of the denominator of specific volume in the Wright EOS [Pa-1]
320
321		! al0 = (a0 + a1T) + a2S
322	0	p0 = (b0 + b4S) + T (b1 + T((b2 + b3T)) + b5*S)
323	0	lambda = (c0 +c4S) + T (c1 + T((c2 + c3T)) + c5*S)
324
325		! SV = al0 + lambda / (pressure + p0)
326
327	0	I_denom = 1.0 / (pressure + p0)
328		dSV_dT = (a1 + I_denom * (c1 + T((2.0c2 + 3.0c3T)) + c5*S)) - &
329	0	(I_denom*2 lambda) * (b1 + T((2.0b2 + 3.0b3T)) + b5*S)
330		dSV_dS = (a2 + I_denom * (c4 + c5*T)) - &
331	0	(I_denom*2 lambda) * (b4 + b5*T)
332
333	0	end subroutine calculate_specvol_derivs_elem_buggy_Wright
334
335		!> Compute the in situ density of sea water (rho) and the compressibility (drho/dp == C_sound^-2)
336		!! at the given salinity, potential temperature and pressure
337		!! using the poor implementation of the equation of state, as fit by Wright, 1997
338	0	elemental subroutine calculate_compress_elem_buggy_Wright(this, T, S, pressure, rho, drho_dp)
339		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
340		real, intent(in) :: T !< Potential temperature relative to the surface [degC]
341		real, intent(in) :: S !< Salinity [PSU]
342		real, intent(in) :: pressure !< Pressure [Pa]
343		real, intent(out) :: rho !< In situ density [kg m-3]
344		real, intent(out) :: drho_dp !< The partial derivative of density with pressure
345		!! (also the inverse of the square of sound speed)
346		!! [s2 m-2].
347
348		! Local variables
349		real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
350		real :: p0 ! The pressure offset in the Wright EOS [Pa]
351		real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
352		real :: I_denom ! The inverse of the denominator of density in the Wright EOS [s2 m-2]
353
354	0	al0 = (a0 + a1T) +a2S
355	0	p0 = (b0 + b4S) + T (b1 + T((b2 + b3T)) + b5*S)
356	0	lambda = (c0 +c4S) + T (c1 + T((c2 + c3T)) + c5*S)
357
358	0	I_denom = 1.0 / (lambda + al0*(pressure + p0))
359	0	rho = (pressure + p0) * I_denom
360	0	drho_dp = lambda * I_denom * I_denom
361
362	0	end subroutine calculate_compress_elem_buggy_Wright
363
364		!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine
365		!! the layer-average specific volumes. There are essentially no free assumptions, apart from a
366		!! truncation in the series for log(1-eps/1+eps) that assumes that \|eps\| < 0.34.
367	0	subroutine avg_spec_vol_buggy_Wright(T, S, p_t, dp, SpV_avg, start, npts)
368		real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface
369		!! [degC].
370		real, dimension(:), intent(in) :: S !< Salinity [PSU].
371		real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa]
372		real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa]
373		real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume
374		!! in the layer [m3 kg-1]
375		integer, intent(in) :: start !< the starting point in the arrays.
376		integer, intent(in) :: npts !< the number of values to calculate.
377
378		! Local variables
379		real :: al0 ! A term in the Wright EOS [m3 kg-1]
380		real :: p0 ! A term in the Wright EOS [Pa]
381		real :: lambda ! A term in the Wright EOS [m2 s-2]
382		real :: eps2 ! The square of a nondimensional ratio [nondim]
383		real :: I_pterm ! The inverse of p0 plus p_ave [Pa-1].
384		real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0, C1_9 = 1.0/9.0 ! Rational constants [nondim]
385		integer :: j
386
387		! alpha(j) = al0 + lambda / (pressure(j) + p0)
388	0	do j=start,start+npts-1
389	0	al0 = a0 + (a1T(j) + a2S(j))
390	0	p0 = b0 + ( b4S(j) + T(j) (b1 + (T(j)(b2 + b3T(j)) + b5*S(j))) )
391	0	lambda = c0 + ( c4S(j) + T(j) (c1 + (T(j)(c2 + c3T(j)) + c5*S(j))) )
392
393	0	I_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j)))
394	0	eps2 = (0.5 * dp(j) * I_pterm)**2
395		SpV_avg(j) = al0 + (lambda * I_pterm) * &
396	0	(1.0 + eps2(C1_3 + eps2(0.2 + eps2(C1_7 + eps2C1_9))))
397		enddo
398	0	end subroutine avg_spec_vol_buggy_Wright
399
400		!> Return the range of temperatures, salinities and pressures for which the reduced-range equation
401		!! of state from Wright (1997) has been fitted to observations. Care should be taken when applying
402		!! this equation of state outside of its fit range.
403	0	subroutine EoS_fit_range_buggy_Wright(this, T_min, T_max, S_min, S_max, p_min, p_max)
404		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
405		real, optional, intent(out) :: T_min !< The minimum potential temperature over which this EoS is fitted [degC]
406		real, optional, intent(out) :: T_max !< The maximum potential temperature over which this EoS is fitted [degC]
407		real, optional, intent(out) :: S_min !< The minimum practical salinity over which this EoS is fitted [PSU]
408		real, optional, intent(out) :: S_max !< The maximum practical salinity over which this EoS is fitted [PSU]
409		real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
410		real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
411
412	0	if (present(T_min)) T_min = -2.0
413	0	if (present(T_max)) T_max = 30.0
414	0	if (present(S_min)) S_min = 28.0
415	0	if (present(S_max)) S_max = 38.0
416	0	if (present(p_min)) p_min = 0.0
417	0	if (present(p_max)) p_max = 5.0e7
418
419	0	end subroutine EoS_fit_range_buggy_Wright
420
421		!> Calculates analytical and nearly-analytical integrals, in geopotential across layers, of pressure
422		!! anomalies, which are required for calculating the finite-volume form pressure accelerations in a
423		!! Boussinesq model. There are essentially no free assumptions, apart from the use of Boole's rule
424		!! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps)
425		!! that assumes that \|eps\| < 0.34.
426	176	subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
427	440	dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, &
428	88	MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p)
429		type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays.
430		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
431		intent(in) :: T !< Potential temperature relative to the surface
432		!! [C ~> degC].
433		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
434		intent(in) :: S !< Salinity [S ~> PSU].
435		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
436		intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m].
437		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
438		intent(in) :: z_b !< Height at the top of the layer [Z ~> m].
439		real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted
440		!! out to reduce the magnitude of each of the integrals.
441		!! (The pressure is calculated as p~=-zrho_0G_e.)
442		real, intent(in) :: rho_0 !< Density [R ~> kg m-3], that is used
443		!! to calculate the pressure (as p~=-zrho_0G_e)
444		!! used in the equation of state.
445		real, intent(in) :: G_e !< The Earth's gravitational acceleration
446		!! [L2 Z-1 T-2 ~> m s-2].
447		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
448		intent(inout) :: dpa !< The change in the pressure anomaly across the
449		!! layer [R L2 T-2 ~> Pa].
450		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
451		optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer
452		!! of the pressure anomaly relative to the anomaly
453		!! at the top of the layer [R Z L2 T-2 ~> Pa m].
454		real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
455		optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
456		!! pressure anomaly at the top and bottom of the
457		!! layer divided by the x grid spacing [R L2 T-2 ~> Pa].
458		real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
459		optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
460		!! pressure anomaly at the top and bottom of the
461		!! layer divided by the y grid spacing [R L2 T-2 ~> Pa].
462		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
463		optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m].
464		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
465		optional, intent(in) :: SSH !< The sea surface height [Z ~> m]
466		real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m].
467		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
468		!! mass weighting to interpolate T/S in integrals
469		real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density
470		!! from kg m-3 to the desired units [R m3 kg-1 ~> 1]
471		real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
472		!! into Pa [Pa T2 R-1 L-2 ~> 1].
473		real, optional, intent(in) :: temp_scale !< A multiplicative factor by which to scale
474		!! temperature into degC [degC C-1 ~> 1]
475		real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure
476		!! into PSU [PSU S-1 ~> 1].
477		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
478		optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
479
480		! Local variables
481	176	real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [m3 kg-1]
482	176	real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: p0_2d ! A term in the Wright EOS [Pa]
483	176	real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: lambda_2d ! A term in the Wright EOS [m2 s-2]
484	176	real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: z0pres ! The height at which the pressure is zero [Z ~> m]
485		real :: al0 ! A term in the Wright EOS [m3 kg-1]
486		real :: p0 ! A term in the Wright EOS [Pa]
487		real :: lambda ! A term in the Wright EOS [m2 s-2]
488		real :: rho_anom ! The density anomaly from rho_ref [kg m-3].
489		real :: eps, eps2 ! A nondimensional ratio and its square [nondim]
490		real :: rem ! [kg m-1 s-2]
491		real :: GxRho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
492		real :: g_Earth ! The gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
493		real :: I_Rho ! The inverse of the Boussinesq density [m3 kg-1]
494		real :: rho_ref_mks ! The reference density in MKS units [kg m-3]
495		real :: p_ave ! The layer averaged pressure [Pa]
496		real :: I_al0 ! The inverse of al0 [kg m-3]
497		real :: I_Lzz ! The inverse of the denominator [Pa-1]
498		real :: dz ! The layer thickness [Z ~> m].
499		real :: hWght ! A pressure-thickness below topography [Z ~> m].
500		real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m].
501		real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2].
502		real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim].
503		real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim].
504		real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim].
505		real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim].
506		real :: intz(5) ! The gravitational acceleration times the integrals of density
507		! with height at the 5 sub-column locations [R L2 T-2 ~> Pa].
508		real :: Pa_to_RL2_T2 ! A conversion factor of pressures from Pa to the output units indicated by
509		! pres_scale [R L2 T-2 Pa-1 ~> 1].
510		real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1]
511		real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1]
512		real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1]
513		real :: b2s ! Partly rescaled version of b2 [Pa C-2 ~> Pa degC-2]
514		real :: b3s ! Partly rescaled version of b3 [Pa C-3 ~> Pa degC-3]
515		real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa PSU-1]
516		real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 PSU-1]
517		real :: c1s ! Partly rescaled version of c1 [m2 s-2 C-1 ~> m2 s-2 degC-1]
518		real :: c2s ! Partly rescaled version of c2 [m2 s-2 C-2 ~> m2 s-2 degC-2]
519		real :: c3s ! Partly rescaled version of c3 [m2 s-2 C-3 ~> m2 s-2 degC-3]
520		real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1]
521		real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1]
522		logical :: do_massWeight ! Indicates whether to do mass weighting.
523		logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
524		real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim]
525		real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim]
526		integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m
527
528		! These array bounds work for the indexing convention of the input arrays, but
529		! on the computational domain defined for the output arrays.
530	88	Isq = HI%IscB ; Ieq = HI%IecB
531	88	Jsq = HI%JscB ; Jeq = HI%JecB
532	88	is = HI%isc ; ie = HI%iec
533	88	js = HI%jsc ; je = HI%jec
534
535	88	if (present(pres_scale)) then
536	0	GxRho = pres_scale * G_e * rho_0 ; g_Earth = pres_scale * G_e
537	0	Pa_to_RL2_T2 = 1.0 / pres_scale
538		else
539	88	GxRho = G_e * rho_0 ; g_Earth = G_e
540	88	Pa_to_RL2_T2 = 1.0
541		endif
542	88	if (present(rho_scale)) then
543	0	g_Earth = g_Earth * rho_scale
544	0	rho_ref_mks = rho_ref / rho_scale ; I_Rho = rho_scale / rho_0
545		else
546	88	rho_ref_mks = rho_ref ; I_Rho = 1.0 / rho_0
547		endif
548		!$omp target enter data map(alloc: z0pres, al0_2d, p0_2d, lambda_2d, intz)
549	88	if (present(Z_0p)) then
550	31944	do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
551	5829736	z0pres(i,j) = Z_0p(i,j)
552		enddo
553		else
554	0	do concurrent (j=HI%jsd:HI%jed, i=HI%isd:HI%ied)
555	0	z0pres(i,j) = 0.0
556		enddo
557		endif
558
559	88	a1s = a1 ; a2s = a2
560	88	b1s = b1 ; b2s = b2 ; b3s = b3 ; b4s = b4 ; b5s = b5
561	88	c1s = c1 ; c2s = c2 ; c3s = c3 ; c4s = c4 ; c5s = c5
562
563	88	if (present(temp_scale)) then ; if (temp_scale /= 1.0) then
564	0	a1s = a1s * temp_scale
565	0	b1s = b1s * temp_scale ; b2s = b2s * temp_scale**2
566	0	b3s = b3s * temp_scale*3 ; b5s = b5s temp_scale
567	0	c1s = c1s * temp_scale ; c2s = c2s * temp_scale**2
568	0	c3s = c3s * temp_scale*3 ; c5s = c5s temp_scale
569		endif ; endif
570
571	88	if (present(saln_scale)) then ; if (saln_scale /= 1.0) then
572	0	a2s = a2s * saln_scale
573	0	b4s = b4s * saln_scale ; b5s = b5s * saln_scale
574	0	c4s = c4s * saln_scale ; c5s = c5s * saln_scale
575		endif ; endif
576
577	88	do_massWeight = .false. ; top_massWeight = .false.
578	88	if (present(MassWghtInterp)) then
579	88	do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
580	88	top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
581		endif
582
583	31944	do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
584	5797792	al0_2d(i,j) = (a0 + a1sT(i,j)) + a2sS(i,j)
585	5797792	p0_2d(i,j) = (b0 + b4sS(i,j)) + T(i,j) (b1s + T(i,j)((b2s + b3sT(i,j))) + b5s*S(i,j))
586	5797792	lambda_2d(i,j) = (c0 +c4sS(i,j)) + T(i,j) (c1s + T(i,j)((c2s + c3sT(i,j))) + c5s*S(i,j))
587
588	5797792	al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
589
590	5797792	dz = z_t(i,j) - z_b(i,j)
591	5797792	p_ave = -GxRho(0.5(z_t(i,j)+z_b(i,j)) - z0pres(i,j))
592
593	5797792	I_al0 = 1.0 / al0
594	5797792	I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)
595	5797792	eps = 0.5GxRhodzI_Lzz ; eps2 = epseps
596
597		! rho(j) = (pressure(j) + p0) / (lambda + al0*(pressure(j) + p0))
598
599	5797792	rho_anom = (p0 + p_ave)(I_LzzI_al0) - rho_ref_mks
600		rem = I_Rho * (lambda * I_al0*2) eps2 * &
601	5797792	(C1_3 + eps2(0.2 + eps2(C1_7 + C1_9*eps2)))
602	5797792	dpa(i,j) = Pa_to_RL2_T2 * (g_Earthrho_anomdz - 2.0epsrem)
603	5797792	if (present(intz_dpa)) &
604	5829736	intz_dpa(i,j) = Pa_to_RL2_T2 * (0.5g_Earthrho_anomdz2 - dz(1.0+eps)*rem)
605		enddo
606
607	88	if (present(intx_dpa)) then
608		!$omp target teams loop collapse(2) &
609		!$omp private(hWght, hL, hR, iDenom, hWt_LL, hWt_LR, hWt_RR, hWt_RL, m, wt_L, wt_R, wtT_L, &
610		!$omp wtT_R, al0, p0, lambda, dz, p_ave, I_al0, I_Lzz, eps, eps2, intz)
611	5734168	do j=js,je ; do I=Isq,Ieq
612		! hWght is the distance measure by which the cell is violation of
613		! hydrostatic consistency. For large hWght we bias the interpolation of
614		! T & S along the top and bottom integrals, akin to thickness weighting.
615	5718240	hWght = 0.0
616	5718240	if (do_massWeight) &
617	0	hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j))
618	5718240	if (top_massWeight) &
619	0	hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j))
620	5718240	if (hWght > 0.) then
621	0	hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
622	0	hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
623	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
624	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
625	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
626	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
627		else
628	5718240	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
629		endif
630
631	5718240	intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
632	22872960	do m=2,4
633	17154720	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
634	17154720	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
635
636	17154720	al0 = (wtT_Lal0_2d(i,j)) + (wtT_Ral0_2d(i+1,j))
637	17154720	p0 = (wtT_Lp0_2d(i,j)) + (wtT_Rp0_2d(i+1,j))
638	17154720	lambda = (wtT_Llambda_2d(i,j)) + (wtT_Rlambda_2d(i+1,j))
639
640	17154720	dz = (wt_L(z_t(i,j) - z_b(i,j))) + (wt_R(z_t(i+1,j) - z_b(i+1,j)))
641		p_ave = -GxRho((wt_L (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + &
642	17154720	(wt_R * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j))))
643
644	17154720	I_al0 = 1.0 / al0
645	17154720	I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)
646	17154720	eps = 0.5GxRhodzI_Lzz ; eps2 = epseps
647
648		intz(m) = Pa_to_RL2_T2 * ( g_Earthdz((p0 + p_ave)(I_LzzI_al0) - rho_ref_mks) - 2.0eps &
649	22872960	I_Rho * (lambda * I_al0*2) eps2 * (C1_3 + eps2(0.2 + eps2(C1_7 + C1_9*eps2))) )
650		enddo
651		! Use Boole's rule to integrate the values.
652	5734080	intx_dpa(i,j) = C1_90(7.0(intz(1)+intz(5)) + 32.0(intz(2)+intz(4)) + 12.0intz(3))
653		enddo ; enddo
654		endif
655
656	88	if (present(inty_dpa)) then
657		!$omp target teams loop collapse(2) &
658		!$omp private(hWght, hL, hR, iDenom, hWt_LL, hWt_LR, hWt_RR, hWt_RL, m, wt_L, wt_R, wtT_L, &
659		!$omp wtT_R, al0, p0, lambda, dz, p_ave, I_al0, I_Lzz, eps, eps2, intz)
660	5750096	do J=Jsq,Jeq ; do i=is,ie
661		! hWght is the distance measure by which the cell is violation of
662		! hydrostatic consistency. For large hWght we bias the interpolation of
663		! T & S along the top and bottom integrals, akin to thickness weighting.
664	5734080	hWght = 0.0
665	5734080	if (do_massWeight) &
666	0	hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j))
667	5734080	if (top_massWeight) &
668	0	hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1))
669	5734080	if (hWght > 0.) then
670	0	hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
671	0	hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
672	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
673	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
674	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
675	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
676		else
677	5734080	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
678		endif
679
680	5734080	intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
681	22936320	do m=2,4
682	17202240	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
683	17202240	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
684
685	17202240	al0 = (wtT_Lal0_2d(i,j)) + (wtT_Ral0_2d(i,j+1))
686	17202240	p0 = (wtT_Lp0_2d(i,j)) + (wtT_Rp0_2d(i,j+1))
687	17202240	lambda = (wtT_Llambda_2d(i,j)) + (wtT_Rlambda_2d(i,j+1))
688
689	17202240	dz = (wt_L(z_t(i,j) - z_b(i,j))) + (wt_R(z_t(i,j+1) - z_b(i,j+1)))
690		p_ave = -GxRho((wt_L(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + &
691	17202240	(wt_R(0.5(z_t(i,j+1)+z_b(i,j+1))-z0pres(i,j+1))))
692
693	17202240	I_al0 = 1.0 / al0
694	17202240	I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)
695	17202240	eps = 0.5GxRhodzI_Lzz ; eps2 = epseps
696
697		intz(m) = Pa_to_RL2_T2 * ( g_Earthdz((p0 + p_ave)(I_LzzI_al0) - rho_ref_mks) - 2.0eps &
698	22936320	I_Rho * (lambda * I_al0*2) eps2 * (C1_3 + eps2(0.2 + eps2(C1_7 + C1_9*eps2))) )
699		enddo
700		! Use Boole's rule to integrate the values.
701	5750008	inty_dpa(i,j) = C1_90(7.0(intz(1)+intz(5)) + 32.0(intz(2)+intz(4)) + 12.0intz(3))
702		enddo ; enddo
703		endif
704		!$omp target exit data map(release: z0pres, al0_2d, p0_2d, lambda_2d, intz)
705
706	528	end subroutine int_density_dz_wright
707
708		!> Calculates analytical and nearly-analytical integrals, in pressure across layers, of geopotential
709		!! anomalies, which are required for calculating the finite-volume form pressure accelerations in a
710		!! non-Boussinesq model. There are essentially no free assumptions, apart from the use of Boole's
711		!! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps)
712		!! that assumes that \|eps\| < 0.34.
713	0	subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
714	0	intp_dza, intx_dza, inty_dza, halo_size, bathyP, P_surf, dP_neglect, &
715		MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale)
716		type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type.
717		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
718		intent(in) :: T !< Potential temperature relative to the surface
719		!! [C ~> degC].
720		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
721		intent(in) :: S !< Salinity [S ~> PSU].
722		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
723		intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
724		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
725		intent(in) :: p_b !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
726		real, intent(in) :: spv_ref !< A mean specific volume that is subtracted out
727		!! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1].
728		!! The calculation is mathematically identical with different values of
729		!! spv_ref, but this reduces the effects of roundoff.
730		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
731		intent(inout) :: dza !< The change in the geopotential anomaly across
732		!! the layer [L2 T-2 ~> m2 s-2].
733		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
734		optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
735		!! the geopotential anomaly relative to the anomaly
736		!! at the bottom of the layer [R L4 T-4 ~> Pa m2 s-2]
737		real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
738		optional, intent(inout) :: intx_dza !< The integral in x of the difference between the
739		!! geopotential anomaly at the top and bottom of
740		!! the layer divided by the x grid spacing
741		!! [L2 T-2 ~> m2 s-2].
742		real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
743		optional, intent(inout) :: inty_dza !< The integral in y of the difference between the
744		!! geopotential anomaly at the top and bottom of
745		!! the layer divided by the y grid spacing
746		!! [L2 T-2 ~> m2 s-2].
747		integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate
748		!! dza.
749		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
750		optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
751		real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
752		optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
753		real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with
754		!! the same units as p_t [R L2 T-2 ~> Pa]
755		integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
756		!! mass weighting to interpolate T/S in integrals
757		real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific
758		!! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
759		real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
760		!! into Pa [Pa T2 R-1 L-2 ~> 1].
761		real, optional, intent(in) :: temp_scale !< A multiplicative factor by which to scale
762		!! temperature into degC [degC C-1 ~> 1]
763		real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure
764		!! into PSU [PSU S-1 ~> 1].
765
766		! Local variables
767	0	real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [R-1 ~> m3 kg-1]
768	0	real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: p0_2d ! A term in the Wright EOS [R L2 T-2 ~> Pa]
769	0	real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: lambda_2d ! A term in the Wright EOS [L2 T-2 ~> m2 s-2]
770		real :: al0 ! A term in the Wright EOS [R-1 ~> m3 kg-1]
771		real :: p0 ! A term in the Wright EOS [R L2 T-2 ~> Pa]
772		real :: lambda ! A term in the Wright EOS [L2 T-2 ~> m2 s-2]
773		real :: al0_scale ! Scaling factor to convert al0 from MKS units [kg m-3 R-1 ~> 1]
774		real :: p0_scale ! Scaling factor to convert p0 from MKS units [R L2 T-2 Pa-1 ~> 1]
775		real :: lam_scale ! Scaling factor to convert lambda from MKS units [L2 s2 T-2 m-2 ~> 1]
776		real :: p_ave ! The layer average pressure [R L2 T-2 ~> Pa]
777		real :: rem ! [L2 T-2 ~> m2 s-2]
778		real :: eps, eps2 ! A nondimensional ratio and its square [nondim]
779		real :: alpha_anom ! The depth averaged specific volume anomaly [R-1 ~> m3 kg-1].
780		real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa].
781		real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa].
782		real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa].
783		real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2].
784		real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim].
785		real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim].
786		real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim].
787		real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim].
788		real :: intp(5) ! The integrals of specific volume with pressure at the
789		! 5 sub-column locations [L2 T-2 ~> m2 s-2].
790		real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1]
791		real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1]
792		real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1]
793		real :: b2s ! Partly rescaled version of b2 [Pa C-2 ~> Pa degC-2]
794		real :: b3s ! Partly rescaled version of b3 [Pa C-3 ~> Pa degC-3]
795		real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa PSU-1]
796		real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 PSU-1]
797		real :: c1s ! Partly rescaled version of c1 [m2 s-2 C-1 ~> m2 s-2 degC-1]
798		real :: c2s ! Partly rescaled version of c2 [m2 s-2 C-2 ~> m2 s-2 degC-2]
799		real :: c3s ! Partly rescaled version of c3 [m2 s-2 C-3 ~> m2 s-2 degC-3]
800		real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1]
801		real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1]
802		logical :: do_massWeight ! Indicates whether to do mass weighting.
803		logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
804		logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting
805		real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim]
806		real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim]
807		integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, halo
808
809	0	Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB
810	0	halo = 0 ; if (present(halo_size)) halo = MAX(halo_size,0)
811	0	ish = HI%isc-halo ; ieh = HI%iec+halo ; jsh = HI%jsc-halo ; jeh = HI%jec+halo
812	0	if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh) ; endif
813	0	if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh) ; endif
814
815
816	0	al0_scale = 1.0 ; if (present(SV_scale)) al0_scale = SV_scale
817	0	p0_scale = 1.0
818	0	if (present(pres_scale)) then ; if (pres_scale /= 1.0) then
819	0	p0_scale = 1.0 / pres_scale
820		endif ; endif
821	0	lam_scale = al0_scale * p0_scale
822
823	0	a1s = a1 ; a2s = a2
824	0	b1s = b1 ; b2s = b2 ; b3s = b3 ; b4s = b4 ; b5s = b5
825	0	c1s = c1 ; c2s = c2 ; c3s = c3 ; c4s = c4 ; c5s = c5
826
827	0	if (present(temp_scale)) then ; if (temp_scale /= 1.0) then
828	0	a1s = a1s * temp_scale
829	0	b1s = b1s * temp_scale ; b2s = b2s * temp_scale**2
830	0	b3s = b3s * temp_scale*3 ; b5s = b5s temp_scale
831	0	c1s = c1s * temp_scale ; c2s = c2s * temp_scale**2
832	0	c3s = c3s * temp_scale*3 ; c5s = c5s temp_scale
833		endif ; endif
834
835	0	if (present(saln_scale)) then ; if (saln_scale /= 1.0) then
836	0	a2s = a2s * saln_scale
837	0	b4s = b4s * saln_scale ; b5s = b5s * saln_scale
838	0	c4s = c4s * saln_scale ; c5s = c5s * saln_scale
839		endif ; endif
840
841	0	do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false.
842	0	if (present(MassWghtInterp)) then
843	0	do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
844	0	top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
845	0	massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set
846		endif
847
848		! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0)
849	0	do j=jsh,jeh ; do i=ish,ieh
850	0	al0_2d(i,j) = al0_scale * ( (a0 + a1sT(i,j)) + a2sS(i,j) )
851	0	p0_2d(i,j) = p0_scale * ( (b0 + b4sS(i,j)) + T(i,j) (b1s + T(i,j)((b2s + b3sT(i,j))) + b5s*S(i,j)) )
852	0	lambda_2d(i,j) = lam_scale * ( (c0 + c4sS(i,j)) + T(i,j) (c1s + T(i,j)((c2s + c3sT(i,j))) + c5s*S(i,j)) )
853
854	0	al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
855	0	dp = p_b(i,j) - p_t(i,j)
856	0	p_ave = 0.5*(p_t(i,j)+p_b(i,j))
857
858	0	eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
859	0	alpha_anom = al0 + lambda / (p0 + p_ave) - spv_ref
860	0	rem = lambda * eps2 * (C1_3 + eps2(0.2 + eps2(C1_7 + C1_9*eps2)))
861	0	dza(i,j) = alpha_anomdp + 2.0eps*rem
862	0	if (present(intp_dza)) &
863	0	intp_dza(i,j) = 0.5alpha_anomdp*2 - dp(1.0-eps)*rem
864		enddo ; enddo
865
866	0	if (present(intx_dza)) then ; do j=HI%jsc,HI%jec ; do I=Isq,Ieq
867		! hWght is the distance measure by which the cell is violation of
868		! hydrostatic consistency. For large hWght we bias the interpolation of
869		! T & S along the top and bottom integrals, akin to thickness weighting.
870	0	hWght = 0.0
871	0	if (do_massWeight .and. massWeight_bug) then
872	0	hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j))
873	0	elseif (do_massWeight) then
874	0	hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j))
875		endif
876	0	if (top_massWeight) &
877	0	hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j))
878	0	if (hWght > 0.) then
879	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
880	0	hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
881	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
882	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
883	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
884	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
885		else
886	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
887		endif
888
889	0	intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
890	0	do m=2,4
891	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
892	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
893
894		! T, S and p are interpolated in the horizontal. The p interpolation
895		! is linear, but for T and S it may be thickness weighted.
896	0	al0 = (wtT_Lal0_2d(i,j)) + (wtT_Ral0_2d(i+1,j))
897	0	p0 = (wtT_Lp0_2d(i,j)) + (wtT_Rp0_2d(i+1,j))
898	0	lambda = (wtT_Llambda_2d(i,j)) + (wtT_Rlambda_2d(i+1,j))
899
900	0	dp = (wt_L(p_b(i,j) - p_t(i,j))) + (wt_R(p_b(i+1,j) - p_t(i+1,j)))
901	0	p_ave = 0.5((wt_L(p_t(i,j)+p_b(i,j))) + (wt_R*(p_t(i+1,j)+p_b(i+1,j))))
902
903	0	eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
904		intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)dp + 2.0eps* &
905	0	lambda * eps2 * (C1_3 + eps2(0.2 + eps2(C1_7 + C1_9*eps2)))
906		enddo
907		! Use Boole's rule to integrate the values.
908		intx_dza(i,j) = C1_90(7.0(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
909	0	12.0*intp(3))
910		enddo ; enddo ; endif
911
912	0	if (present(inty_dza)) then ; do J=Jsq,Jeq ; do i=HI%isc,HI%iec
913		! hWght is the distance measure by which the cell is violation of
914		! hydrostatic consistency. For large hWght we bias the interpolation of
915		! T & S along the top and bottom integrals, akin to thickness weighting.
916	0	hWght = 0.0
917	0	if (do_massWeight .and. massWeight_bug) then
918	0	hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j))
919	0	elseif (do_massWeight) then
920	0	hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1))
921		endif
922	0	if (top_massWeight) &
923	0	hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j))
924	0	if (hWght > 0.) then
925	0	hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
926	0	hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
927	0	hWght = hWght * ( (hL-hR)/(hL+hR) )**2
928	0	iDenom = 1.0 / ( hWght(hR + hL) + hLhR )
929	0	hWt_LL = (hWghthL + hRhL) * iDenom ; hWt_LR = (hWghthR) iDenom
930	0	hWt_RR = (hWghthR + hRhL) * iDenom ; hWt_RL = (hWghthL) iDenom
931		else
932	0	hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
933		endif
934
935	0	intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
936	0	do m=2,4
937	0	wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
938	0	wtT_L = (wt_LhWt_LL) + (wt_RhWt_RL) ; wtT_R = (wt_LhWt_LR) + (wt_RhWt_RR)
939
940		! T, S and p are interpolated in the horizontal. The p interpolation
941		! is linear, but for T and S it may be thickness weighted.
942	0	al0 = (wt_Lal0_2d(i,j)) + (wt_Ral0_2d(i,j+1))
943	0	p0 = (wt_Lp0_2d(i,j)) + (wt_Rp0_2d(i,j+1))
944	0	lambda = (wt_Llambda_2d(i,j)) + (wt_Rlambda_2d(i,j+1))
945
946	0	dp = (wt_L(p_b(i,j) - p_t(i,j))) + (wt_R(p_b(i,j+1) - p_t(i,j+1)))
947	0	p_ave = 0.5((wt_L(p_t(i,j)+p_b(i,j))) + (wt_R*(p_t(i,j+1)+p_b(i,j+1))))
948
949	0	eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
950		intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)dp + 2.0eps* &
951	0	lambda * eps2 * (C1_3 + eps2(0.2 + eps2(C1_7 + C1_9*eps2)))
952		enddo
953		! Use Boole's rule to integrate the values.
954		inty_dza(i,j) = C1_90(7.0(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
955	0	12.0*intp(3))
956		enddo ; enddo ; endif
957	0	end subroutine int_spec_vol_dp_wright
958
959		!> Calculate the in-situ density for 1D arraya inputs and outputs.
960	17393	subroutine calculate_density_array_buggy_Wright(this, T, S, pressure, rho, start, npts, rho_ref)
961		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
962		real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
963		real, dimension(:), intent(in) :: S !< Salinity [PSU]
964		real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
965		real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]
966		integer, intent(in) :: start !< The starting index for calculations
967		integer, intent(in) :: npts !< The number of values to calculate
968		real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]
969
970		! Local variables
971		integer :: j
972
973	17393	if (present(rho_ref)) then
974	0	do j = start, start+npts-1
975	0	rho(j) = density_anomaly_elem_buggy_Wright(this, T(j), S(j), pressure(j), rho_ref)
976		enddo
977		else
978	6278843	do j = start, start+npts-1
979	6278843	rho(j) = density_elem_buggy_Wright_loc(T(j), S(j), pressure(j))
980		enddo
981		endif
982	17393	end subroutine calculate_density_array_buggy_Wright
983
984		!> Calculate the in-situ density for 2D arraya inputs and outputs.
985	12	subroutine calculate_density_array_2d_buggy_Wright(this, T, S, pressure, rho, &
986		dom, rho_ref)
987		class(buggy_Wright_EOS), intent(in) :: this
988		!< This EOS
989		real, intent(in) :: T(:,:)
990		!< Potential temperature relative to the surface [degC]
991		real, intent(in) :: S(:,:)
992		!< Salinity [PSU]
993		real, intent(in) :: pressure(:,:)
994		!< Pressure [Pa]
995		real, intent(out) :: rho(:,:)
996		!< In situ density [kg m-3]
997		integer, intent(in) :: dom(2,2)
998		!< Index bounds of domain. First index is rank, second is bounds
999		real, optional, intent(in) :: rho_ref
1000		!< A reference density [kg m-3]
1001
1002		integer :: is, ie, js, je
1003		integer :: i, j
1004
1005	12	is = dom(1,1) ; ie = dom(1,2)
1006	12	js = dom(2,1) ; je = dom(2,2)
1007
1008		! NOTE: There is an implicit copy of `this` which cannot yet be prevented.
1009		! Possibly because Nvidia cannot associate `this` with `EOS%type`.
1010
1011	12	if (present(rho_ref)) then
1012	0	do concurrent (j=js:je, i=is:ie)
1013		rho(i,j) = density_anomaly_elem_buggy_Wright(this, T(i,j), S(i,j), &
1014	0	pressure(i,j), rho_ref)
1015		enddo
1016		else
1017	12	do concurrent (j=js:je, i=is:ie)
1018	808188	rho(i,j) = density_elem_buggy_Wright_loc( T(i,j), S(i,j), pressure(i,j))
1019		enddo
1020		endif
1021	12	end subroutine calculate_density_array_2d_buggy_Wright
1022
1023		!> Calculate the in-situ density for 3D array inputs and outputs.
1024	0	subroutine calculate_density_array_3d_buggy_Wright(this, T, S, pressure, rho, &
1025		dom, rho_ref)
1026		class(buggy_Wright_EOS), intent(in) :: this
1027		!< This EOS
1028		real, intent(in) :: T(:,:,:)
1029		!< Potential temperature relative to the surface [degC]
1030		real, intent(in) :: S(:,:,:)
1031		!< Salinity [PSU]
1032		real, intent(in) :: pressure(:,:,:)
1033		!< Pressure [Pa]
1034		real, intent(out) :: rho(:,:,:)
1035		!< In situ density [kg m-3]
1036		integer, intent(in) :: dom(3,2)
1037		!< Index bounds of domain. First index is rank, second is bounds
1038		real, optional, intent(in) :: rho_ref
1039		!< A reference density [kg m-3]
1040
1041		integer :: is, ie, js, je, ks, ke
1042		integer :: i, j, k
1043
1044	0	is = dom(1,1) ; ie = dom(1,2)
1045	0	js = dom(2,1) ; je = dom(2,2)
1046	0	ks = dom(3,1) ; ke = dom(3,2)
1047
1048		! NOTE: There is an implicit copy of `this` which cannot yet be prevented.
1049		! Possibly because Nvidia cannot associate `this` with `EOS%type`.
1050
1051	0	if (present(rho_ref)) then
1052	0	do concurrent (k=ks:ke, j=js:je, i=is:ie)
1053		rho(i,j,k) = density_anomaly_elem_buggy_Wright(this, T(i,j,k), S(i,j,k), &
1054	0	pressure(i,j,k), rho_ref)
1055		enddo
1056		else
1057	0	do concurrent (k=ks:ke, j=js:je, i=is:ie)
1058	0	rho(i,j,k) = density_elem_buggy_Wright_loc( T(i,j,k), S(i,j,k), pressure(i,j,k))
1059		enddo
1060		endif
1061	0	end subroutine calculate_density_array_3d_buggy_Wright
1062
1063		!> Calculate the in-situ specific volume for 1D array inputs and outputs.
1064	0	subroutine calculate_spec_vol_array_buggy_Wright(this, T, S, pressure, specvol, start, npts, spv_ref)
1065		class(buggy_Wright_EOS), intent(in) :: this !< This EOS
1066		real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
1067		real, dimension(:), intent(in) :: S !< Salinity [PSU]
1068		real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
1069		real, dimension(:), intent(out) :: specvol !< In situ specific volume [m3 kg-1]
1070		integer, intent(in) :: start !< The starting index for calculations
1071		integer, intent(in) :: npts !< The number of values to calculate
1072		real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
1073
1074		! Local variables
1075		integer :: j
1076
1077	0	if (present(spv_ref)) then
1078	0	do j = start, start+npts-1
1079	0	specvol(j) = spec_vol_anomaly_elem_buggy_Wright(this, T(j), S(j), pressure(j), spv_ref)
1080		enddo
1081		else
1082	0	do j = start, start+npts-1
1083	0	specvol(j) = spec_vol_elem_buggy_Wright(this, T(j), S(j), pressure(j) )
1084		enddo
1085		endif
1086
1087	0	end subroutine calculate_spec_vol_array_buggy_Wright
1088
1089
1090		!> Calculate the in-situ density derivatives for 2D array inputs and outputs.
1091	0	subroutine calculate_density_derivs_2d_buggy_Wright(this, T, S, pressure, &
1092	0	drho_dT, drho_dS, dom)
1093		class(buggy_Wright_EOS), intent(in) :: this
1094		!< This EOS
1095		real, intent(in) :: T(:,:)
1096		!< Potential temperature relative to the surface [degC]
1097		real, intent(in) :: S(:,:)
1098		!< Salinity [PSU]
1099		real, intent(in) :: pressure(:,:)
1100		!< Pressure [Pa]
1101		real, intent(out) :: drho_dT(:,:)
1102		!< Partial derivative of density with potential temperature [kg m-3 degC-1]
1103		real, intent(out) :: drho_dS(:,:)
1104		!< Partial derivative of density with salinity [kg m-3 PSU-1]
1105		integer, intent(in) :: dom(2,2)
1106		!< Index bounds of domain. First index is rank, second is bounds
1107
1108		integer :: is, ie, js, je
1109		integer :: i, j
1110
1111	0	is = dom(1,1) ; ie = dom(1,2)
1112	0	js = dom(2,1) ; je = dom(2,2)
1113
1114		! NOTE: There is an implicit copy of `this` which cannot yet be prevented.
1115
1116	0	do concurrent (j=js:je, i=is:ie)
1117		call calculate_density_derivs_elem_buggy_Wright_loc( T(i,j), S(i,j), &
1118	0	pressure(i,j), drho_dT(i,j), drho_dS(i,j))
1119		enddo
1120	0	end subroutine calculate_density_derivs_2d_buggy_Wright
1121
1122		!> Calculate the in-situ density derivatives for 3D array inputs and outputs.
1123	84	subroutine calculate_density_derivs_3d_buggy_Wright(this, T, S, pressure, &
1124	84	drho_dT, drho_dS, dom)
1125		class(buggy_Wright_EOS), intent(in) :: this
1126		!< This EOS
1127		real, intent(in) :: T(:,:,:)
1128		!< Potential temperature relative to the surface [degC]
1129		real, intent(in) :: S(:,:,:)
1130		!< Salinity [PSU]
1131		real, intent(in) :: pressure(:,:,:)
1132		!< Pressure [Pa]
1133		real, intent(out) :: drho_dT(:,:,:)
1134		!< Partial derivative of density with potential temperature [kg m-3 degC-1]
1135		real, intent(out) :: drho_dS(:,:,:)
1136		!< Partial derivative of density with salinity [kg m-3 PSU-1]
1137		integer, intent(in) :: dom(3,2)
1138		!< Index bounds of domain. First index is rank, second is bounds
1139
1140		integer :: is, ie, js, je, ks, ke
1141		integer :: i, j, k
1142
1143	84	is = dom(1,1) ; ie = dom(1,2)
1144	84	js = dom(2,1) ; je = dom(2,2)
1145	84	ks = dom(3,1) ; ke = dom(3,2)
1146
1147		! NOTE: There is an implicit copy of `this` which cannot yet be prevented.
1148
1149	84	do concurrent (k=ks:ke, j=js:je, i=is:ie)
1150		call calculate_density_derivs_elem_buggy_Wright_loc( T(i,j,k), S(i,j,k), &
1151	11099004	pressure(i,j,k), drho_dT(i,j,k), drho_dS(i,j,k))
1152		enddo
1153	84	end subroutine calculate_density_derivs_3d_buggy_Wright
1154
1155		!> Set coefficients that can correct bugs un the buggy Wright equation of state.
1156	2	subroutine set_params_buggy_Wright(this, use_Wright_2nd_deriv_bug)
1157		class(buggy_Wright_EOS), intent(inout) :: this !< This EOS
1158		logical, optional, intent(in) :: use_Wright_2nd_deriv_bug !< If true, use a buggy
1159		!! buggy version of the calculations of the second
1160		!! derivative of density with temperature and with temperature and
1161		!! pressure. This bug is corrected in the default version.
1162
1163	2	this%three = 3.0
1164	2	if (present(use_Wright_2nd_deriv_bug)) then
1165	1	if (use_Wright_2nd_deriv_bug) then ; this%three = 2.0
1166	1	else ; this%three = 3.0 ; endif
1167		endif
1168
1169	2	end subroutine set_params_buggy_Wright
1170
1171
1172		!> \namespace mom_eos_wright
1173		!!
1174		!! \section section_EOS_Wright Wright equation of state
1175		!!
1176		!! Wright, 1997, provide an approximation for the in situ density as a function of
1177		!! potential temperature, salinity and pressure. The formula follow the Tumlirz
1178		!! equation of state which are easier to evaluate and make efficient.
1179		!!
1180		!! Two ranges are provided by Wright: a "full" range and "reduced" range. The version in this
1181		!! module uses the reduced range.
1182		!!
1183		!! Originally coded in 2000 by R. Hallberg.
1184		!! Anomaly form coded in 3/18.
1185		!!
1186		!! \subsection section_EOS_Wright_references References
1187		!!
1188		!! Wright, D., 1997: An Equation of State for Use in Ocean Models: Eckart's Formula Revisited.
1189		!! J. Ocean. Atmosph. Tech., 14 (3), 735-740.
1190		!! https://journals.ametsoc.org/doi/abs/10.1175/1520-0426%281997%29014%3C0735%3AAEOSFU%3E2.0.CO%3B2
1191
1192	2	end module MOM_EOS_Wright