MOM_EOS_Wright.F90

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
8
11
12implicit none ; private
13
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]
25real, parameter :: a0 = 7.057924e-4 ! A parameter in the Wright alpha_0 fit [m3 kg-1]
26real, parameter :: a1 = 3.480336e-7 ! A parameter in the Wright alpha_0 fit [m3 kg-1 degC-1]
27real, parameter :: a2 = -1.112733e-7 ! A parameter in the Wright alpha_0 fit [m3 kg-1 PSU-1]
28real, parameter :: b0 = 5.790749e8 ! A parameter in the Wright p_0 fit [Pa]
29real, parameter :: b1 = 3.516535e6 ! A parameter in the Wright p_0 fit [Pa degC-1]
30real, parameter :: b2 = -4.002714e4 ! A parameter in the Wright p_0 fit [Pa degC-2]
31real, parameter :: b3 = 2.084372e2 ! A parameter in the Wright p_0 fit [Pa degC-3]
32real, parameter :: b4 = 5.944068e5 ! A parameter in the Wright p_0 fit [Pa PSU-1]
33real, parameter :: b5 = -9.643486e3 ! A parameter in the Wright p_0 fit [Pa degC-1 PSU-1]
34real, parameter :: c0 = 1.704853e5 ! A parameter in the Wright lambda fit [m2 s-2]
35real, parameter :: c1 = 7.904722e2 ! A parameter in the Wright lambda fit [m2 s-2 degC-1]
36real, parameter :: c2 = -7.984422 ! A parameter in the Wright lambda fit [m2 s-2 degC-2]
37real, parameter :: c3 = 5.140652e-2 ! A parameter in the Wright lambda fit [m2 s-2 degC-3]
38real, parameter :: c4 = -2.302158e2 ! A parameter in the Wright lambda fit [m2 s-2 PSU-1]
39real, 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
43type, extends (eos_base) :: buggy_wright_eos
44
45 real :: three = 3.0 !< A constant that can be adjusted to recreate some bugs [nondim]
46
47contains
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
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_spec_vol_array for efficiency
73 procedure :: calculate_spec_vol_array => calculate_spec_vol_array_buggy_wright
74
75end type buggy_wright_eos
76
77contains
78
79!> In situ density of sea water using a buggy implementation of Wright, 1997 [kg m-3]
80!!
81!! This is an elemental function that can be applied to any combination of scalar and array inputs.
82real elemental function density_elem_buggy_wright(this, t, s, pressure)
83 class(buggy_wright_eos), intent(in) :: this !< This EOS
84 real, intent(in) :: t !< potential temperature relative to the surface [degC].
85 real, intent(in) :: s !< salinity [PSU].
86 real, intent(in) :: pressure !< pressure [Pa].
87
88 ! Local variables
89 real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
90 real :: p0 ! The pressure offset in the Wright EOS [Pa]
91 real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
92
93 al0 = (a0 + a1*t) +a2*s
94 p0 = (b0 + b4*s) + t * (b1 + t*(b2 + b3*t) + b5*s)
95 lambda = (c0 +c4*s) + t * (c1 + t*(c2 + c3*t) + c5*s)
96 density_elem_buggy_wright = (pressure + p0) / (lambda + al0*(pressure + p0))
97
99
100!> In situ density anomaly of sea water using a buggy implementation of Wright, 1997 [kg m-3]
101!!
102!! This is an elemental function that can be applied to any combination of scalar and array inputs.
103real elemental function density_anomaly_elem_buggy_wright(this, t, s, pressure, rho_ref)
104 class(buggy_wright_eos), intent(in) :: this !< This EOS
105 real, intent(in) :: t !< potential temperature relative to the surface [degC].
106 real, intent(in) :: s !< salinity [PSU].
107 real, intent(in) :: pressure !< pressure [Pa].
108 real, intent(in) :: rho_ref !< A reference density [kg m-3].
109
110 ! Local variables
111 real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
112 real :: p0 ! The pressure offset in the Wright EOS [Pa]
113 real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
114 real :: al_ts ! The contributions of temperature and salinity to al0 [m3 kg-1]
115 real :: p_tsp ! A combination of the pressure and the temperature and salinity contributions to p0 [Pa]
116 real :: lam_ts ! The contributions of temperature and salinity to lambda [m2 s-2]
117 real :: pa_000 ! A corrected offset to the pressure, including contributions from rho_ref [Pa]
118
119 pa_000 = (b0*(1.0 - a0*rho_ref) - rho_ref*c0)
120 al_ts = a1*t +a2*s
121 al0 = a0 + al_ts
122 p_tsp = pressure + (b4*s + t * (b1 + (t*(b2 + b3*t) + b5*s)))
123 lam_ts = c4*s + t * (c1 + (t*(c2 + c3*t) + c5*s))
124
125 ! The following two expressions are mathematically equivalent.
126 ! wright_density = (b0 + p0_TSp) / ((c0 + lam_TS) + al0*(b0 + p0_TSp)) - rho_ref
128 (pa_000 + (p_tsp - rho_ref*(p_tsp*al0 + (b0*al_ts + lam_ts)))) / &
129 ( (c0 + lam_ts) + al0*(b0 + p_tsp) )
130
132
133!> In situ specific volume of sea water using a buggy implementation of Wright, 1997 [kg m-3]
134!!
135!! This is an elemental function that can be applied to any combination of scalar and array inputs.
136real elemental function spec_vol_elem_buggy_wright(this, t, s, pressure)
137 class(buggy_wright_eos), intent(in) :: this !< This EOS
138 real, intent(in) :: t !< potential temperature relative to the surface [degC].
139 real, intent(in) :: s !< salinity [PSU].
140 real, intent(in) :: pressure !< pressure [Pa].
141
142 ! Local variables
143 real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
144 real :: p0 ! The pressure offset in the Wright EOS [Pa]
145 real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
146
147 al0 = (a0 + a1*t) +a2*s
148 p0 = (b0 + b4*s) + t * (b1 + t*((b2 + b3*t)) + b5*s)
149 lambda = (c0 +c4*s) + t * (c1 + t*((c2 + c3*t)) + c5*s)
150
151 spec_vol_elem_buggy_wright = (lambda + al0*(pressure + p0)) / (pressure + p0)
152
154
155!> In situ specific volume anomaly of sea water using a buggy implementation of Wright, 1997 [kg m-3]
156!!
157!! This is an elemental function that can be applied to any combination of scalar and array inputs.
158real elemental function spec_vol_anomaly_elem_buggy_wright(this, t, s, pressure, spv_ref)
159 class(buggy_wright_eos), intent(in) :: this !< This EOS
160 real, intent(in) :: t !< potential temperature relative to the surface [degC].
161 real, intent(in) :: s !< salinity [PSU].
162 real, intent(in) :: pressure !< pressure [Pa].
163 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
164
165 ! Local variables
166 real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
167 real :: p0 ! The pressure offset in the Wright EOS [Pa]
168 real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
169
170 al0 = (a0 + a1*t) +a2*s
171 p0 = (b0 + b4*s) + t * (b1 + t*((b2 + b3*t)) + b5*s)
172 lambda = (c0 +c4*s) + t * (c1 + t*((c2 + c3*t)) + c5*s)
173
174 spec_vol_anomaly_elem_buggy_wright = (lambda + (al0 - spv_ref)*(pressure + p0)) / (pressure + p0)
175
177
178!> Calculate the partial derivatives of density with potential temperature and salinity
179!! using the buggy implementation of the equation of state, as fit by Wright, 1997
180elemental subroutine calculate_density_derivs_elem_buggy_wright(this, T, S, pressure, drho_dT, drho_dS)
181 class(buggy_wright_eos), intent(in) :: this !< This EOS
182 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
183 real, intent(in) :: s !< Salinity [PSU]
184 real, intent(in) :: pressure !< Pressure [Pa]
185 real, intent(out) :: drho_dt !< The partial derivative of density with potential
186 !! temperature [kg m-3 degC-1]
187 real, intent(out) :: drho_ds !< The partial derivative of density with salinity,
188 !! in [kg m-3 PSU-1]
189 ! Local variables
190 real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
191 real :: p0 ! The pressure offset in the Wright EOS [Pa]
192 real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
193 real :: i_denom2 ! The inverse of the square of the denominator of density in the Wright EOS [s4 m-4]
194
195 al0 = (a0 + a1*t) + a2*s
196 p0 = (b0 + b4*s) + t * (b1 + t*((b2 + b3*t)) + b5*s)
197 lambda = (c0 +c4*s) + t * (c1 + t*((c2 + c3*t)) + c5*s)
198
199 i_denom2 = 1.0 / (lambda + al0*(pressure + p0))
200 i_denom2 = i_denom2 *i_denom2
201 drho_dt = i_denom2 * &
202 (lambda* (b1 + t*(2.0*b2 + 3.0*b3*t) + b5*s) - &
203 (pressure+p0) * ( (pressure+p0)*a1 + &
204 (c1 + t*(c2*2.0 + c3*3.0*t) + c5*s) ))
205 drho_ds = i_denom2 * (lambda* (b4 + b5*t) - &
206 (pressure+p0) * ( (pressure+p0)*a2 + (c4 + c5*t) ))
207
209
210!> Second derivatives of density with respect to temperature, salinity, and pressure,
211!! using the poor implementation of the equation of state, as fit by Wright, 1997
212elemental subroutine calculate_density_second_derivs_elem_buggy_wright(this, T, S, pressure, &
213 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
214 class(buggy_wright_eos), intent(in) :: this !< This EOS
215 real, intent(in) :: t !< Potential temperature referenced to 0 dbar [degC]
216 real, intent(in) :: s !< Salinity [PSU]
217 real, intent(in) :: pressure !< Pressure [Pa]
218 real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
219 !! to S [kg m-3 PSU-2]
220 real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
221 !! to T [kg m-3 PSU-1 degC-1]
222 real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
223 !! to T [kg m-3 degC-2]
224 real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
225 !! to pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1]
226 real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
227 !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
228 ! Local variables
229 real :: z0, z1 ! Local work variables [Pa]
230 real :: z2, z4 ! Local work variables [m2 s-2]
231 real :: z3, z5 ! Local work variables [Pa degC-1]
232 real :: z6, z8 ! Local work variables [m2 s-2 degC-1]
233 real :: z7 ! A local work variable [m2 s-2 PSU-1]
234 real :: z9 ! A local work variable [m3 kg-1]
235 real :: z10 ! A local work variable [Pa PSU-1]
236 real :: z11 ! A local work variable [Pa m2 s-2 PSU-1] = [kg m s-4 PSU-1]
237 real :: z2_2 ! A local work variable [m4 s-4]
238 real :: z2_3 ! A local work variable [m6 s-6]
239 real :: six ! A constant that can be adjusted from 6. to 4. to recreate a bug [nondim]
240
241 ! Based on the above expression with common terms factored, there probably exists a more numerically stable
242 ! and/or efficient expression
243
244 six = 2.0*this%three ! When recreating a bug from the original version of this routine, six = 4.
245
246 z0 = t*(b1 + b5*s + t*(b2 + b3*t))
247 z1 = (b0 + pressure + b4*s + z0)
248 z3 = (b1 + b5*s + t*(2.*b2 + this%three*b3*t)) ! When recreating a bug here this%three = 2.
249 z4 = (c0 + c4*s + t*(c1 + c5*s + t*(c2 + c3*t)))
250 z5 = (b1 + b5*s + t*(b2 + b3*t) + t*(b2 + 2.*b3*t))
251 z6 = c1 + c5*s + t*(c2 + c3*t) + t*(c2 + 2.*c3*t)
252 z7 = (c4 + c5*t + a2*z1)
253 z8 = (c1 + c5*s + t*(2.*c2 + 3.*c3*t) + a1*z1)
254 z9 = (a0 + a2*s + a1*t)
255 z10 = (b4 + b5*t)
256 z11 = (z10*z4 - z1*z7)
257 z2 = (c0 + c4*s + t*(c1 + c5*s + t*(c2 + c3*t)) + z9*z1)
258 z2_2 = z2*z2
259 z2_3 = z2_2*z2
260
261 drho_ds_ds = (z10*(c4 + c5*t) - a2*z10*z1 - z10*z7)/z2_2 - (2.*(c4 + c5*t + z9*z10 + a2*z1)*z11)/z2_3
262 drho_ds_dt = (z10*z6 - z1*(c5 + a2*z5) + b5*z4 - z5*z7)/z2_2 - (2.*(z6 + z9*z5 + a1*z1)*z11)/z2_3
263 drho_dt_dt = (z3*z6 - z1*(2.*c2 + 6.*c3*t + a1*z5) + (2.*b2 + six*b3*t)*z4 - z5*z8)/z2_2 - &
264 (2.*(z6 + z9*z5 + a1*z1)*(z3*z4 - z1*z8))/z2_3
265 drho_ds_dp = (-c4 - c5*t - 2.*a2*z1)/z2_2 - (2.*z9*z11)/z2_3
266 drho_dt_dp = (-c1 - c5*s - t*(2.*c2 + 3.*c3*t) - 2.*a1*z1)/z2_2 - (2.*z9*(z3*z4 - z1*z8))/z2_3
267
269
270!> Calculate the partial derivatives of specific volume with temperature and salinity
271!! using the poor implementation of the equation of state, as fit by Wright, 1997
272elemental subroutine calculate_specvol_derivs_elem_buggy_wright(this, T, S, pressure, dSV_dT, dSV_dS)
273 class(buggy_wright_eos), intent(in) :: this !< This EOS
274 real, intent(in) :: t !< Potential temperature [degC]
275 real, intent(in) :: s !< Salinity [PSU]
276 real, intent(in) :: pressure !< Pressure [Pa]
277 real, intent(inout) :: dsv_dt !< The partial derivative of specific volume with
278 !! potential temperature [m3 kg-1 degC-1]
279 real, intent(inout) :: dsv_ds !< The partial derivative of specific volume with
280 !! salinity [m3 kg-1 PSU-1]
281 ! Local variables
282 real :: p0 ! The pressure offset in the Wright EOS [Pa]
283 real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
284 real :: i_denom ! The inverse of the denominator of specific volume in the Wright EOS [Pa-1]
285
286! al0 = (a0 + a1*T) + a2*S
287 p0 = (b0 + b4*s) + t * (b1 + t*((b2 + b3*t)) + b5*s)
288 lambda = (c0 +c4*s) + t * (c1 + t*((c2 + c3*t)) + c5*s)
289
290 ! SV = al0 + lambda / (pressure + p0)
291
292 i_denom = 1.0 / (pressure + p0)
293 dsv_dt = (a1 + i_denom * (c1 + t*((2.0*c2 + 3.0*c3*t)) + c5*s)) - &
294 (i_denom**2 * lambda) * (b1 + t*((2.0*b2 + 3.0*b3*t)) + b5*s)
295 dsv_ds = (a2 + i_denom * (c4 + c5*t)) - &
296 (i_denom**2 * lambda) * (b4 + b5*t)
297
299
300!> Compute the in situ density of sea water (rho) and the compressibility (drho/dp == C_sound^-2)
301!! at the given salinity, potential temperature and pressure
302!! using the poor implementation of the equation of state, as fit by Wright, 1997
303elemental subroutine calculate_compress_elem_buggy_wright(this, T, S, pressure, rho, drho_dp)
304 class(buggy_wright_eos), intent(in) :: this !< This EOS
305 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
306 real, intent(in) :: s !< Salinity [PSU]
307 real, intent(in) :: pressure !< Pressure [Pa]
308 real, intent(out) :: rho !< In situ density [kg m-3]
309 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
310 !! (also the inverse of the square of sound speed)
311 !! [s2 m-2].
312
313 ! Local variables
314 real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1]
315 real :: p0 ! The pressure offset in the Wright EOS [Pa]
316 real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2]
317 real :: i_denom ! The inverse of the denominator of density in the Wright EOS [s2 m-2]
318
319 al0 = (a0 + a1*t) +a2*s
320 p0 = (b0 + b4*s) + t * (b1 + t*((b2 + b3*t)) + b5*s)
321 lambda = (c0 +c4*s) + t * (c1 + t*((c2 + c3*t)) + c5*s)
322
323 i_denom = 1.0 / (lambda + al0*(pressure + p0))
324 rho = (pressure + p0) * i_denom
325 drho_dp = lambda * i_denom * i_denom
326
328
329!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine
330!! the layer-average specific volumes. There are essentially no free assumptions, apart from a
331!! truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
332subroutine avg_spec_vol_buggy_wright(T, S, p_t, dp, SpV_avg, start, npts)
333 real, dimension(:), intent(in) :: t !< Potential temperature relative to the surface
334 !! [degC].
335 real, dimension(:), intent(in) :: s !< Salinity [PSU].
336 real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa]
337 real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa]
338 real, dimension(:), intent(inout) :: spv_avg !< The vertical average specific volume
339 !! in the layer [m3 kg-1]
340 integer, intent(in) :: start !< the starting point in the arrays.
341 integer, intent(in) :: npts !< the number of values to calculate.
342
343 ! Local variables
344 real :: al0 ! A term in the Wright EOS [m3 kg-1]
345 real :: p0 ! A term in the Wright EOS [Pa]
346 real :: lambda ! A term in the Wright EOS [m2 s-2]
347 real :: eps2 ! The square of a nondimensional ratio [nondim]
348 real :: i_pterm ! The inverse of p0 plus p_ave [Pa-1].
349 real, parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0, c1_9 = 1.0/9.0 ! Rational constants [nondim]
350 integer :: j
351
352 ! alpha(j) = al0 + lambda / (pressure(j) + p0)
353 do j=start,start+npts-1
354 al0 = a0 + (a1*t(j) + a2*s(j))
355 p0 = b0 + ( b4*s(j) + t(j) * (b1 + (t(j)*(b2 + b3*t(j)) + b5*s(j))) )
356 lambda = c0 + ( c4*s(j) + t(j) * (c1 + (t(j)*(c2 + c3*t(j)) + c5*s(j))) )
357
358 i_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j)))
359 eps2 = (0.5 * dp(j) * i_pterm)**2
360 spv_avg(j) = al0 + (lambda * i_pterm) * &
361 (1.0 + eps2*(c1_3 + eps2*(0.2 + eps2*(c1_7 + eps2*c1_9))))
362 enddo
363end subroutine avg_spec_vol_buggy_wright
364
365!> Return the range of temperatures, salinities and pressures for which the reduced-range equation
366!! of state from Wright (1997) has been fitted to observations. Care should be taken when applying
367!! this equation of state outside of its fit range.
368subroutine eos_fit_range_buggy_wright(this, T_min, T_max, S_min, S_max, p_min, p_max)
369 class(buggy_wright_eos), intent(in) :: this !< This EOS
370 real, optional, intent(out) :: T_min !< The minimum potential temperature over which this EoS is fitted [degC]
371 real, optional, intent(out) :: T_max !< The maximum potential temperature over which this EoS is fitted [degC]
372 real, optional, intent(out) :: S_min !< The minimum practical salinity over which this EoS is fitted [PSU]
373 real, optional, intent(out) :: S_max !< The maximum practical salinity over which this EoS is fitted [PSU]
374 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
375 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
376
377 if (present(t_min)) t_min = -2.0
378 if (present(t_max)) t_max = 30.0
379 if (present(s_min)) s_min = 28.0
380 if (present(s_max)) s_max = 38.0
381 if (present(p_min)) p_min = 0.0
382 if (present(p_max)) p_max = 5.0e7
383
384end subroutine eos_fit_range_buggy_wright
385
386!> Calculates analytical and nearly-analytical integrals, in geopotential across layers, of pressure
387!! anomalies, which are required for calculating the finite-volume form pressure accelerations in a
388!! Boussinesq model. There are essentially no free assumptions, apart from the use of Boole's rule
389!! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps)
390!! that assumes that |eps| < 0.34.
391subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
392 dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, &
393 MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p)
394 type(hor_index_type), intent(in) :: hi !< The horizontal index type for the arrays.
395 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
396 intent(in) :: t !< Potential temperature relative to the surface
397 !! [C ~> degC].
398 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
399 intent(in) :: s !< Salinity [S ~> PSU].
400 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
401 intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m].
402 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
403 intent(in) :: z_b !< Height at the top of the layer [Z ~> m].
404 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted
405 !! out to reduce the magnitude of each of the integrals.
406 !! (The pressure is calculated as p~=-z*rho_0*G_e.)
407 real, intent(in) :: rho_0 !< Density [R ~> kg m-3], that is used
408 !! to calculate the pressure (as p~=-z*rho_0*G_e)
409 !! used in the equation of state.
410 real, intent(in) :: g_e !< The Earth's gravitational acceleration
411 !! [L2 Z-1 T-2 ~> m s-2].
412 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
413 intent(inout) :: dpa !< The change in the pressure anomaly across the
414 !! layer [R L2 T-2 ~> Pa].
415 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
416 optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer
417 !! of the pressure anomaly relative to the anomaly
418 !! at the top of the layer [R Z L2 T-2 ~> Pa m].
419 real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
420 optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
421 !! pressure anomaly at the top and bottom of the
422 !! layer divided by the x grid spacing [R L2 T-2 ~> Pa].
423 real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
424 optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
425 !! pressure anomaly at the top and bottom of the
426 !! layer divided by the y grid spacing [R L2 T-2 ~> Pa].
427 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
428 optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m].
429 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
430 optional, intent(in) :: ssh !< The sea surface height [Z ~> m]
431 real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m].
432 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use
433 !! mass weighting to interpolate T/S in integrals
434 real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density
435 !! from kg m-3 to the desired units [R m3 kg-1 ~> 1]
436 real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
437 !! into Pa [Pa T2 R-1 L-2 ~> 1].
438 real, optional, intent(in) :: temp_scale !< A multiplicative factor by which to scale
439 !! temperature into degC [degC C-1 ~> 1]
440 real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure
441 !! into PSU [PSU S-1 ~> 1].
442 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
443 optional, intent(in) :: z_0p !< The height at which the pressure is 0 [Z ~> m]
444
445 ! Local variables
446 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [m3 kg-1]
447 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: p0_2d ! A term in the Wright EOS [Pa]
448 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: lambda_2d ! A term in the Wright EOS [m2 s-2]
449 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: z0pres ! The height at which the pressure is zero [Z ~> m]
450 real :: al0 ! A term in the Wright EOS [m3 kg-1]
451 real :: p0 ! A term in the Wright EOS [Pa]
452 real :: lambda ! A term in the Wright EOS [m2 s-2]
453 real :: rho_anom ! The density anomaly from rho_ref [kg m-3].
454 real :: eps, eps2 ! A nondimensional ratio and its square [nondim]
455 real :: rem ! [kg m-1 s-2]
456 real :: gxrho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
457 real :: g_earth ! The gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
458 real :: i_rho ! The inverse of the Boussinesq density [m3 kg-1]
459 real :: rho_ref_mks ! The reference density in MKS units [kg m-3]
460 real :: p_ave ! The layer averaged pressure [Pa]
461 real :: i_al0 ! The inverse of al0 [kg m-3]
462 real :: i_lzz ! The inverse of the denominator [Pa-1]
463 real :: dz ! The layer thickness [Z ~> m].
464 real :: hwght ! A pressure-thickness below topography [Z ~> m].
465 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m].
466 real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2].
467 real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
468 real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
469 real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
470 real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
471 real :: intz(5) ! The gravitational acceleration times the integrals of density
472 ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa].
473 real :: pa_to_rl2_t2 ! A conversion factor of pressures from Pa to the output units indicated by
474 ! pres_scale [R L2 T-2 Pa-1 ~> 1].
475 real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1]
476 real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1]
477 real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1]
478 real :: b2s ! Partly rescaled version of b2 [Pa C-2 ~> Pa degC-2]
479 real :: b3s ! Partly rescaled version of b3 [Pa C-3 ~> Pa degC-3]
480 real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa PSU-1]
481 real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 PSU-1]
482 real :: c1s ! Partly rescaled version of c1 [m2 s-2 C-1 ~> m2 s-2 degC-1]
483 real :: c2s ! Partly rescaled version of c2 [m2 s-2 C-2 ~> m2 s-2 degC-2]
484 real :: c3s ! Partly rescaled version of c3 [m2 s-2 C-3 ~> m2 s-2 degC-3]
485 real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1]
486 real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1]
487 logical :: do_massweight ! Indicates whether to do mass weighting.
488 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface
489 real, parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0 ! Rational constants [nondim]
490 real, parameter :: c1_9 = 1.0/9.0, c1_90 = 1.0/90.0 ! Rational constants [nondim]
491 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m
492
493 ! These array bounds work for the indexing convention of the input arrays, but
494 ! on the computational domain defined for the output arrays.
495 isq = hi%IscB ; ieq = hi%IecB
496 jsq = hi%JscB ; jeq = hi%JecB
497 is = hi%isc ; ie = hi%iec
498 js = hi%jsc ; je = hi%jec
499
500 if (present(pres_scale)) then
501 gxrho = pres_scale * g_e * rho_0 ; g_earth = pres_scale * g_e
502 pa_to_rl2_t2 = 1.0 / pres_scale
503 else
504 gxrho = g_e * rho_0 ; g_earth = g_e
505 pa_to_rl2_t2 = 1.0
506 endif
507 if (present(rho_scale)) then
508 g_earth = g_earth * rho_scale
509 rho_ref_mks = rho_ref / rho_scale ; i_rho = rho_scale / rho_0
510 else
511 rho_ref_mks = rho_ref ; i_rho = 1.0 / rho_0
512 endif
513 if (present(z_0p)) then
514 do j=jsq,jeq+1 ; do i=isq,ieq+1
515 z0pres(i,j) = z_0p(i,j)
516 enddo ; enddo
517 else
518 z0pres(:,:) = 0.0
519 endif
520
521 a1s = a1 ; a2s = a2
522 b1s = b1 ; b2s = b2 ; b3s = b3 ; b4s = b4 ; b5s = b5
523 c1s = c1 ; c2s = c2 ; c3s = c3 ; c4s = c4 ; c5s = c5
524
525 if (present(temp_scale)) then ; if (temp_scale /= 1.0) then
526 a1s = a1s * temp_scale
527 b1s = b1s * temp_scale ; b2s = b2s * temp_scale**2
528 b3s = b3s * temp_scale**3 ; b5s = b5s * temp_scale
529 c1s = c1s * temp_scale ; c2s = c2s * temp_scale**2
530 c3s = c3s * temp_scale**3 ; c5s = c5s * temp_scale
531 endif ; endif
532
533 if (present(saln_scale)) then ; if (saln_scale /= 1.0) then
534 a2s = a2s * saln_scale
535 b4s = b4s * saln_scale ; b5s = b5s * saln_scale
536 c4s = c4s * saln_scale ; c5s = c5s * saln_scale
537 endif ; endif
538
539 do_massweight = .false. ; top_massweight = .false.
540 if (present(masswghtinterp)) then
541 do_massweight = btest(masswghtinterp, 0) ! True for odd values
542 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set
543 endif
544
545 do j=jsq,jeq+1 ; do i=isq,ieq+1
546 al0_2d(i,j) = (a0 + a1s*t(i,j)) + a2s*s(i,j)
547 p0_2d(i,j) = (b0 + b4s*s(i,j)) + t(i,j) * (b1s + t(i,j)*((b2s + b3s*t(i,j))) + b5s*s(i,j))
548 lambda_2d(i,j) = (c0 +c4s*s(i,j)) + t(i,j) * (c1s + t(i,j)*((c2s + c3s*t(i,j))) + c5s*s(i,j))
549
550 al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
551
552 dz = z_t(i,j) - z_b(i,j)
553 p_ave = -gxrho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))
554
555 i_al0 = 1.0 / al0
556 i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
557 eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
558
559! rho(j) = (pressure(j) + p0) / (lambda + al0*(pressure(j) + p0))
560
561 rho_anom = (p0 + p_ave)*(i_lzz*i_al0) - rho_ref_mks
562 rem = i_rho * (lambda * i_al0**2) * eps2 * &
563 (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
564 dpa(i,j) = pa_to_rl2_t2 * (g_earth*rho_anom*dz - 2.0*eps*rem)
565 if (present(intz_dpa)) &
566 intz_dpa(i,j) = pa_to_rl2_t2 * (0.5*g_earth*rho_anom*dz**2 - dz*(1.0+eps)*rem)
567 enddo ; enddo
568
569 if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
570 ! hWght is the distance measure by which the cell is violation of
571 ! hydrostatic consistency. For large hWght we bias the interpolation of
572 ! T & S along the top and bottom integrals, akin to thickness weighting.
573 hwght = 0.0
574 if (do_massweight) &
575 hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
576 if (top_massweight) &
577 hwght = max(hwght, z_b(i+1,j)-ssh(i,j), z_b(i,j)-ssh(i+1,j))
578 if (hwght > 0.) then
579 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
580 hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
581 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
582 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
583 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
584 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
585 else
586 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
587 endif
588
589 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
590 do m=2,4
591 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
592 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
593
594 al0 = (wtt_l*al0_2d(i,j)) + (wtt_r*al0_2d(i+1,j))
595 p0 = (wtt_l*p0_2d(i,j)) + (wtt_r*p0_2d(i+1,j))
596 lambda = (wtt_l*lambda_2d(i,j)) + (wtt_r*lambda_2d(i+1,j))
597
598 dz = (wt_l*(z_t(i,j) - z_b(i,j))) + (wt_r*(z_t(i+1,j) - z_b(i+1,j)))
599 p_ave = -gxrho*((wt_l * (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + &
600 (wt_r * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j))))
601
602 i_al0 = 1.0 / al0
603 i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
604 eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
605
606 intz(m) = pa_to_rl2_t2 * ( g_earth*dz*((p0 + p_ave)*(i_lzz*i_al0) - rho_ref_mks) - 2.0*eps * &
607 i_rho * (lambda * i_al0**2) * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2))) )
608 enddo
609 ! Use Boole's rule to integrate the values.
610 intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))
611 enddo ; enddo ; endif
612
613 if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
614 ! hWght is the distance measure by which the cell is violation of
615 ! hydrostatic consistency. For large hWght we bias the interpolation of
616 ! T & S along the top and bottom integrals, akin to thickness weighting.
617 hwght = 0.0
618 if (do_massweight) &
619 hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
620 if (top_massweight) &
621 hwght = max(hwght, z_b(i,j+1)-ssh(i,j), z_b(i,j)-ssh(i,j+1))
622 if (hwght > 0.) then
623 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
624 hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
625 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
626 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
627 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
628 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
629 else
630 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
631 endif
632
633 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
634 do m=2,4
635 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
636 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
637
638 al0 = (wtt_l*al0_2d(i,j)) + (wtt_r*al0_2d(i,j+1))
639 p0 = (wtt_l*p0_2d(i,j)) + (wtt_r*p0_2d(i,j+1))
640 lambda = (wtt_l*lambda_2d(i,j)) + (wtt_r*lambda_2d(i,j+1))
641
642 dz = (wt_l*(z_t(i,j) - z_b(i,j))) + (wt_r*(z_t(i,j+1) - z_b(i,j+1)))
643 p_ave = -gxrho*((wt_l*(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + &
644 (wt_r*(0.5*(z_t(i,j+1)+z_b(i,j+1))-z0pres(i,j+1))))
645
646 i_al0 = 1.0 / al0
647 i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
648 eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
649
650 intz(m) = pa_to_rl2_t2 * ( g_earth*dz*((p0 + p_ave)*(i_lzz*i_al0) - rho_ref_mks) - 2.0*eps * &
651 i_rho * (lambda * i_al0**2) * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2))) )
652 enddo
653 ! Use Boole's rule to integrate the values.
654 inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))
655 enddo ; enddo ; endif
656
657end subroutine int_density_dz_wright
658
659!> Calculates analytical and nearly-analytical integrals, in pressure across layers, of geopotential
660!! anomalies, which are required for calculating the finite-volume form pressure accelerations in a
661!! non-Boussinesq model. There are essentially no free assumptions, apart from the use of Boole's
662!! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps)
663!! that assumes that |eps| < 0.34.
664subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
665 intp_dza, intx_dza, inty_dza, halo_size, bathyP, P_surf, dP_neglect, &
666 MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale)
667 type(hor_index_type), intent(in) :: hi !< The ocean's horizontal index type.
668 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
669 intent(in) :: t !< Potential temperature relative to the surface
670 !! [C ~> degC].
671 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
672 intent(in) :: s !< Salinity [S ~> PSU].
673 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
674 intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
675 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
676 intent(in) :: p_b !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
677 real, intent(in) :: spv_ref !< A mean specific volume that is subtracted out
678 !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1].
679 !! The calculation is mathematically identical with different values of
680 !! spv_ref, but this reduces the effects of roundoff.
681 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
682 intent(inout) :: dza !< The change in the geopotential anomaly across
683 !! the layer [L2 T-2 ~> m2 s-2].
684 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
685 optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
686 !! the geopotential anomaly relative to the anomaly
687 !! at the bottom of the layer [R L4 T-4 ~> Pa m2 s-2]
688 real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
689 optional, intent(inout) :: intx_dza !< The integral in x of the difference between the
690 !! geopotential anomaly at the top and bottom of
691 !! the layer divided by the x grid spacing
692 !! [L2 T-2 ~> m2 s-2].
693 real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
694 optional, intent(inout) :: inty_dza !< The integral in y of the difference between the
695 !! geopotential anomaly at the top and bottom of
696 !! the layer divided by the y grid spacing
697 !! [L2 T-2 ~> m2 s-2].
698 integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate
699 !! dza.
700 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
701 optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
702 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
703 optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
704 real, optional, intent(in) :: dp_neglect !< A miniscule pressure change with
705 !! the same units as p_t [R L2 T-2 ~> Pa]
706 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use
707 !! mass weighting to interpolate T/S in integrals
708 real, optional, intent(in) :: sv_scale !< A multiplicative factor by which to scale specific
709 !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
710 real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
711 !! into Pa [Pa T2 R-1 L-2 ~> 1].
712 real, optional, intent(in) :: temp_scale !< A multiplicative factor by which to scale
713 !! temperature into degC [degC C-1 ~> 1]
714 real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure
715 !! into PSU [PSU S-1 ~> 1].
716
717 ! Local variables
718 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [R-1 ~> m3 kg-1]
719 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: p0_2d ! A term in the Wright EOS [R L2 T-2 ~> Pa]
720 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: lambda_2d ! A term in the Wright EOS [L2 T-2 ~> m2 s-2]
721 real :: al0 ! A term in the Wright EOS [R-1 ~> m3 kg-1]
722 real :: p0 ! A term in the Wright EOS [R L2 T-2 ~> Pa]
723 real :: lambda ! A term in the Wright EOS [L2 T-2 ~> m2 s-2]
724 real :: al0_scale ! Scaling factor to convert al0 from MKS units [kg m-3 R-1 ~> 1]
725 real :: p0_scale ! Scaling factor to convert p0 from MKS units [R L2 T-2 Pa-1 ~> 1]
726 real :: lam_scale ! Scaling factor to convert lambda from MKS units [L2 s2 T-2 m-2 ~> 1]
727 real :: p_ave ! The layer average pressure [R L2 T-2 ~> Pa]
728 real :: rem ! [L2 T-2 ~> m2 s-2]
729 real :: eps, eps2 ! A nondimensional ratio and its square [nondim]
730 real :: alpha_anom ! The depth averaged specific volume anomaly [R-1 ~> m3 kg-1].
731 real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa].
732 real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa].
733 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa].
734 real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2].
735 real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
736 real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
737 real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
738 real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
739 real :: intp(5) ! The integrals of specific volume with pressure at the
740 ! 5 sub-column locations [L2 T-2 ~> m2 s-2].
741 real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1]
742 real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1]
743 real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1]
744 real :: b2s ! Partly rescaled version of b2 [Pa C-2 ~> Pa degC-2]
745 real :: b3s ! Partly rescaled version of b3 [Pa C-3 ~> Pa degC-3]
746 real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa PSU-1]
747 real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 PSU-1]
748 real :: c1s ! Partly rescaled version of c1 [m2 s-2 C-1 ~> m2 s-2 degC-1]
749 real :: c2s ! Partly rescaled version of c2 [m2 s-2 C-2 ~> m2 s-2 degC-2]
750 real :: c3s ! Partly rescaled version of c3 [m2 s-2 C-3 ~> m2 s-2 degC-3]
751 real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1]
752 real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1]
753 logical :: do_massweight ! Indicates whether to do mass weighting.
754 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface
755 logical :: massweight_bug ! If true, use an incorrect expression to determine where to apply mass weighting
756 real, parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0 ! Rational constants [nondim]
757 real, parameter :: c1_9 = 1.0/9.0, c1_90 = 1.0/90.0 ! Rational constants [nondim]
758 integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, halo
759
760 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
761 halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
762 ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
763 if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh) ; endif
764 if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh) ; endif
765
766
767 al0_scale = 1.0 ; if (present(sv_scale)) al0_scale = sv_scale
768 p0_scale = 1.0
769 if (present(pres_scale)) then ; if (pres_scale /= 1.0) then
770 p0_scale = 1.0 / pres_scale
771 endif ; endif
772 lam_scale = al0_scale * p0_scale
773
774 a1s = a1 ; a2s = a2
775 b1s = b1 ; b2s = b2 ; b3s = b3 ; b4s = b4 ; b5s = b5
776 c1s = c1 ; c2s = c2 ; c3s = c3 ; c4s = c4 ; c5s = c5
777
778 if (present(temp_scale)) then ; if (temp_scale /= 1.0) then
779 a1s = a1s * temp_scale
780 b1s = b1s * temp_scale ; b2s = b2s * temp_scale**2
781 b3s = b3s * temp_scale**3 ; b5s = b5s * temp_scale
782 c1s = c1s * temp_scale ; c2s = c2s * temp_scale**2
783 c3s = c3s * temp_scale**3 ; c5s = c5s * temp_scale
784 endif ; endif
785
786 if (present(saln_scale)) then ; if (saln_scale /= 1.0) then
787 a2s = a2s * saln_scale
788 b4s = b4s * saln_scale ; b5s = b5s * saln_scale
789 c4s = c4s * saln_scale ; c5s = c5s * saln_scale
790 endif ; endif
791
792 do_massweight = .false. ; massweight_bug = .false. ; top_massweight = .false.
793 if (present(masswghtinterp)) then
794 do_massweight = btest(masswghtinterp, 0) ! True for odd values
795 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set
796 massweight_bug = btest(masswghtinterp, 3) ! True if the 8 bit is set
797 endif
798
799 ! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0)
800 do j=jsh,jeh ; do i=ish,ieh
801 al0_2d(i,j) = al0_scale * ( (a0 + a1s*t(i,j)) + a2s*s(i,j) )
802 p0_2d(i,j) = p0_scale * ( (b0 + b4s*s(i,j)) + t(i,j) * (b1s + t(i,j)*((b2s + b3s*t(i,j))) + b5s*s(i,j)) )
803 lambda_2d(i,j) = lam_scale * ( (c0 + c4s*s(i,j)) + t(i,j) * (c1s + t(i,j)*((c2s + c3s*t(i,j))) + c5s*s(i,j)) )
804
805 al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
806 dp = p_b(i,j) - p_t(i,j)
807 p_ave = 0.5*(p_t(i,j)+p_b(i,j))
808
809 eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
810 alpha_anom = al0 + lambda / (p0 + p_ave) - spv_ref
811 rem = lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
812 dza(i,j) = alpha_anom*dp + 2.0*eps*rem
813 if (present(intp_dza)) &
814 intp_dza(i,j) = 0.5*alpha_anom*dp**2 - dp*(1.0-eps)*rem
815 enddo ; enddo
816
817 if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
818 ! hWght is the distance measure by which the cell is violation of
819 ! hydrostatic consistency. For large hWght we bias the interpolation of
820 ! T & S along the top and bottom integrals, akin to thickness weighting.
821 hwght = 0.0
822 if (do_massweight .and. massweight_bug) then
823 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
824 elseif (do_massweight) then
825 hwght = max(0., p_t(i+1,j)-bathyp(i,j), p_t(i,j)-bathyp(i+1,j))
826 endif
827 if (top_massweight) &
828 hwght = max(hwght, p_surf(i,j)-p_b(i+1,j), p_surf(i+1,j)-p_b(i,j))
829 if (hwght > 0.) then
830 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
831 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
832 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
833 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
834 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
835 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
836 else
837 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
838 endif
839
840 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
841 do m=2,4
842 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
843 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
844
845 ! T, S and p are interpolated in the horizontal. The p interpolation
846 ! is linear, but for T and S it may be thickness weighted.
847 al0 = (wtt_l*al0_2d(i,j)) + (wtt_r*al0_2d(i+1,j))
848 p0 = (wtt_l*p0_2d(i,j)) + (wtt_r*p0_2d(i+1,j))
849 lambda = (wtt_l*lambda_2d(i,j)) + (wtt_r*lambda_2d(i+1,j))
850
851 dp = (wt_l*(p_b(i,j) - p_t(i,j))) + (wt_r*(p_b(i+1,j) - p_t(i+1,j)))
852 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))))
853
854 eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
855 intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
856 lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
857 enddo
858 ! Use Boole's rule to integrate the values.
859 intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
860 12.0*intp(3))
861 enddo ; enddo ; endif
862
863 if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
864 ! hWght is the distance measure by which the cell is violation of
865 ! hydrostatic consistency. For large hWght we bias the interpolation of
866 ! T & S along the top and bottom integrals, akin to thickness weighting.
867 hwght = 0.0
868 if (do_massweight .and. massweight_bug) then
869 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
870 elseif (do_massweight) then
871 hwght = max(0., p_t(i,j+1)-bathyp(i,j), p_t(i,j)-bathyp(i,j+1))
872 endif
873 if (top_massweight) &
874 hwght = max(hwght, p_surf(i,j)-p_b(i,j+1), p_surf(i,j+1)-p_b(i,j))
875 if (hwght > 0.) then
876 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
877 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
878 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
879 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
880 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
881 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
882 else
883 hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
884 endif
885
886 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
887 do m=2,4
888 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
889 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
890
891 ! T, S and p are interpolated in the horizontal. The p interpolation
892 ! is linear, but for T and S it may be thickness weighted.
893 al0 = (wt_l*al0_2d(i,j)) + (wt_r*al0_2d(i,j+1))
894 p0 = (wt_l*p0_2d(i,j)) + (wt_r*p0_2d(i,j+1))
895 lambda = (wt_l*lambda_2d(i,j)) + (wt_r*lambda_2d(i,j+1))
896
897 dp = (wt_l*(p_b(i,j) - p_t(i,j))) + (wt_r*(p_b(i,j+1) - p_t(i,j+1)))
898 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))))
899
900 eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
901 intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
902 lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
903 enddo
904 ! Use Boole's rule to integrate the values.
905 inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
906 12.0*intp(3))
907 enddo ; enddo ; endif
908end subroutine int_spec_vol_dp_wright
909
910!> Calculate the in-situ density for 1D arraya inputs and outputs.
911subroutine calculate_density_array_buggy_wright(this, T, S, pressure, rho, start, npts, rho_ref)
912 class(buggy_wright_eos), intent(in) :: this !< This EOS
913 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
914 real, dimension(:), intent(in) :: S !< Salinity [PSU]
915 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
916 real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]
917 integer, intent(in) :: start !< The starting index for calculations
918 integer, intent(in) :: npts !< The number of values to calculate
919 real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]
920
921 ! Local variables
922 integer :: j
923
924 if (present(rho_ref)) then
925 do j = start, start+npts-1
926 rho(j) = density_anomaly_elem_buggy_wright(this, t(j), s(j), pressure(j), rho_ref)
927 enddo
928 else
929 do j = start, start+npts-1
930 rho(j) = density_elem_buggy_wright(this, t(j), s(j), pressure(j))
931 enddo
932 endif
933
935
936!> Calculate the in-situ specific volume for 1D array inputs and outputs.
937subroutine calculate_spec_vol_array_buggy_wright(this, T, S, pressure, specvol, start, npts, spv_ref)
938 class(buggy_wright_eos), intent(in) :: this !< This EOS
939 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
940 real, dimension(:), intent(in) :: S !< Salinity [PSU]
941 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
942 real, dimension(:), intent(out) :: specvol !< In situ specific volume [m3 kg-1]
943 integer, intent(in) :: start !< The starting index for calculations
944 integer, intent(in) :: npts !< The number of values to calculate
945 real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
946
947 ! Local variables
948 integer :: j
949
950 if (present(spv_ref)) then
951 do j = start, start+npts-1
952 specvol(j) = spec_vol_anomaly_elem_buggy_wright(this, t(j), s(j), pressure(j), spv_ref)
953 enddo
954 else
955 do j = start, start+npts-1
956 specvol(j) = spec_vol_elem_buggy_wright(this, t(j), s(j), pressure(j) )
957 enddo
958 endif
959
961
962
963!> Set coefficients that can correct bugs un the buggy Wright equation of state.
964subroutine set_params_buggy_wright(this, use_Wright_2nd_deriv_bug)
965 class(buggy_wright_eos), intent(inout) :: this !< This EOS
966 logical, optional, intent(in) :: use_wright_2nd_deriv_bug !< If true, use a buggy
967 !! buggy version of the calculations of the second
968 !! derivative of density with temperature and with temperature and
969 !! pressure. This bug is corrected in the default version.
970
971 this%three = 3.0
972 if (present(use_wright_2nd_deriv_bug)) then
973 if (use_wright_2nd_deriv_bug) then ; this%three = 2.0
974 else ; this%three = 3.0 ; endif
975 endif
976
977end subroutine set_params_buggy_wright
978
979
980!> \namespace mom_eos_wright
981!!
982!! \section section_EOS_Wright Wright equation of state
983!!
984!! Wright, 1997, provide an approximation for the in situ density as a function of
985!! potential temperature, salinity and pressure. The formula follow the Tumlirz
986!! equation of state which are easier to evaluate and make efficient.
987!!
988!! Two ranges are provided by Wright: a "full" range and "reduced" range. The version in this
989!! module uses the reduced range.
990!!
991!! Originally coded in 2000 by R. Hallberg.
992!! Anomaly form coded in 3/18.
993!!
994!! \subsection section_EOS_Wright_references References
995!!
996!! Wright, D., 1997: An Equation of State for Use in Ocean Models: Eckart's Formula Revisited.
997!! J. Ocean. Atmosph. Tech., 14 (3), 735-740.
998!! https://journals.ametsoc.org/doi/abs/10.1175/1520-0426%281997%29014%3C0735%3AAEOSFU%3E2.0.CO%3B2
999
1000end module mom_eos_wright