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