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