MOM_EOS_UNESCO.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 Jackett and McDougall fits to the UNESCO EOS
7
9
10implicit none ; private
11
12public unesco_eos
13
14!>@{ Parameters in the UNESCO equation of state, as published in appendix A3 of Gill, 1982.
15! The following constants are used to calculate rho0, the density of seawater at 1 atmosphere pressure.
16! The notation is Rab for the contribution to rho0 from S^a*T^b, with 6 used for the 1.5 power.
17real, parameter :: r00 = 999.842594 ! A coefficient in the fit for rho0 [kg m-3]
18real, parameter :: r01 = 6.793952e-2 ! A coefficient in the fit for rho0 [kg m-3 degC-1]
19real, parameter :: r02 = -9.095290e-3 ! A coefficient in the fit for rho0 [kg m-3 degC-2]
20real, parameter :: r03 = 1.001685e-4 ! A coefficient in the fit for rho0 [kg m-3 degC-3]
21real, parameter :: r04 = -1.120083e-6 ! A coefficient in the fit for rho0 [kg m-3 degC-4]
22real, parameter :: r05 = 6.536332e-9 ! A coefficient in the fit for rho0 [kg m-3 degC-5]
23real, parameter :: r10 = 0.824493 ! A coefficient in the fit for rho0 [kg m-3 PSU-1]
24real, parameter :: r11 = -4.0899e-3 ! A coefficient in the fit for rho0 [kg m-3 degC-1 PSU-1]
25real, parameter :: r12 = 7.6438e-5 ! A coefficient in the fit for rho0 [kg m-3 degC-2 PSU-1]
26real, parameter :: r13 = -8.2467e-7 ! A coefficient in the fit for rho0 [kg m-3 degC-3 PSU-1]
27real, parameter :: r14 = 5.3875e-9 ! A coefficient in the fit for rho0 [kg m-3 degC-4 PSU-1]
28real, parameter :: r60 = -5.72466e-3 ! A coefficient in the fit for rho0 [kg m-3 PSU-1.5]
29real, parameter :: r61 = 1.0227e-4 ! A coefficient in the fit for rho0 [kg m-3 degC-1 PSU-1.5]
30real, parameter :: r62 = -1.6546e-6 ! A coefficient in the fit for rho0 [kg m-3 degC-2 PSU-1.5]
31real, parameter :: r20 = 4.8314e-4 ! A coefficient in the fit for rho0 [kg m-3 PSU-2]
32
33! The following constants are used to calculate the secant bulk modulus.
34! The notation here is Sabc for terms proportional to S^a*T^b*P^c, with 6 used for the 1.5 power.
35! Note that these values differ from those in Appendix 3 of Gill (1982) because the expressions
36! from Jackett and MacDougall (1995) use potential temperature, rather than in situ temperature.
37real, parameter :: s000 = 1.965933e4 ! A coefficient in the secant bulk modulus fit [bar]
38real, parameter :: s010 = 1.444304e2 ! A coefficient in the secant bulk modulus fit [bar degC-1]
39real, parameter :: s020 = -1.706103 ! A coefficient in the secant bulk modulus fit [bar degC-2]
40real, parameter :: s030 = 9.648704e-3 ! A coefficient in the secant bulk modulus fit [bar degC-3]
41real, parameter :: s040 = -4.190253e-5 ! A coefficient in the secant bulk modulus fit [bar degC-4]
42real, parameter :: s100 = 52.84855 ! A coefficient in the secant bulk modulus fit [bar PSU-1]
43real, parameter :: s110 = -3.101089e-1 ! A coefficient in the secant bulk modulus fit [bar degC-1 PSU-1]
44real, parameter :: s120 = 6.283263e-3 ! A coefficient in the secant bulk modulus fit [bar degC-2 PSU-1]
45real, parameter :: s130 = -5.084188e-5 ! A coefficient in the secant bulk modulus fit [bar degC-3 PSU-1]
46real, parameter :: s600 = 3.886640e-1 ! A coefficient in the secant bulk modulus fit [bar PSU-1.5]
47real, parameter :: s610 = 9.085835e-3 ! A coefficient in the secant bulk modulus fit [bar degC-1 PSU-1.5]
48real, parameter :: s620 = -4.619924e-4 ! A coefficient in the secant bulk modulus fit [bar degC-2 PSU-1.5]
49
50real, parameter :: s001 = 3.186519 ! A coefficient in the secant bulk modulus fit [nondim]
51real, parameter :: s011 = 2.212276e-2 ! A coefficient in the secant bulk modulus fit [degC-1]
52real, parameter :: s021 = -2.984642e-4 ! A coefficient in the secant bulk modulus fit [degC-2]
53real, parameter :: s031 = 1.956415e-6 ! A coefficient in the secant bulk modulus fit [degC-3]
54real, parameter :: s101 = 6.704388e-3 ! A coefficient in the secant bulk modulus fit [PSU-1]
55real, parameter :: s111 = -1.847318e-4 ! A coefficient in the secant bulk modulus fit [degC-1 PSU-1]
56real, parameter :: s121 = 2.059331e-7 ! A coefficient in the secant bulk modulus fit [degC-2 PSU-1]
57real, parameter :: s601 = 1.480266e-4 ! A coefficient in the secant bulk modulus fit [PSU-1.5]
58
59real, parameter :: s002 = 2.102898e-4 ! A coefficient in the secant bulk modulus fit [bar-1]
60real, parameter :: s012 = -1.202016e-5 ! A coefficient in the secant bulk modulus fit [bar-1 degC-1]
61real, parameter :: s022 = 1.394680e-7 ! A coefficient in the secant bulk modulus fit [bar-1 degC-2]
62real, parameter :: s102 = -2.040237e-6 ! A coefficient in the secant bulk modulus fit [bar-1 PSU-1]
63real, parameter :: s112 = 6.128773e-8 ! A coefficient in the secant bulk modulus fit [bar-1 degC-1 PSU-1]
64real, parameter :: s122 = 6.207323e-10 ! A coefficient in the secant bulk modulus fit [bar-1 degC-2 PSU-1]
65!>@}
66
67!> The EOS_base implementation of the UNESCO equation of state
68type, extends (eos_base) :: unesco_eos
69
70contains
71 !> Implementation of the in-situ density as an elemental function [kg m-3]
72 procedure :: density_elem => density_elem_unesco
73 !> Implementation of the in-situ density anomaly as an elemental function [kg m-3]
74 procedure :: density_anomaly_elem => density_anomaly_elem_unesco
75 !> Implementation of the in-situ specific volume as an elemental function [m3 kg-1]
76 procedure :: spec_vol_elem => spec_vol_elem_unesco
77 !> Implementation of the in-situ specific volume anomaly as an elemental function [m3 kg-1]
78 procedure :: spec_vol_anomaly_elem => spec_vol_anomaly_elem_unesco
79 !> Implementation of the calculation of derivatives of density
80 procedure :: calculate_density_derivs_elem => calculate_density_derivs_elem_unesco
81 !> Implementation of the calculation of second derivatives of density
82 procedure :: calculate_density_second_derivs_elem => calculate_density_second_derivs_elem_unesco
83 !> Implementation of the calculation of derivatives of specific volume
84 procedure :: calculate_specvol_derivs_elem => calculate_specvol_derivs_elem_unesco
85 !> Implementation of the calculation of compressibility
86 procedure :: calculate_compress_elem => calculate_compress_elem_unesco
87 !> Implementation of the range query function
88 procedure :: eos_fit_range => eos_fit_range_unesco
89
90end type unesco_eos
91
92contains
93
94!> In situ density as fit by Jackett and McDougall, 1995 [kg m-3]
95!!
96!! This is an elemental function that can be applied to any combination of scalar and array inputs.
97real elemental function density_elem_unesco(this, t, s, pressure)
98 class(unesco_eos), intent(in) :: this !< This EOS
99 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
100 real, intent(in) :: s !< Salinity [PSU]
101 real, intent(in) :: pressure !< Pressure [Pa]
102
103 ! Local variables
104 real :: t1 ! A copy of the temperature at a point [degC]
105 real :: s1 ! A copy of the salinity at a point [PSU]
106 real :: p1 ! Pressure converted to bars [bar]
107 real :: s12 ! The square root of salinity [PSU1/2]
108 real :: rho0 ! Density at 1 bar pressure [kg m-3]
109 real :: sig0 ! The anomaly of rho0 from R00 [kg m-3]
110 real :: ks ! The secant bulk modulus [bar]
111
112 p1 = pressure*1.0e-5 ; t1 = t
113 s1 = max(s, 0.0) ; s12 = sqrt(s1)
114
115 ! Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ).
116 sig0 = ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
117 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
118 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
119 rho0 = r00 + sig0
120
121 ! Compute rho(s,theta,p), first calculating the secant bulk modulus.
122 ks = (s000 + ( t1*(s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
123 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620))) )) + &
124 p1*( (s001 + ( t1*(s011 + t1*(s021 + t1*s031)) + &
125 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )) + &
126 p1*(s002 + ( t1*(s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )) )
127
128 density_elem_unesco = rho0*ks / (ks - p1)
129
130end function density_elem_unesco
131
132!> In situ density anomaly as fit by Jackett and McDougall, 1995 [kg m-3]
133!!
134!! This is an elemental function that can be applied to any combination of scalar and array inputs.
135real elemental function density_anomaly_elem_unesco(this, t, s, pressure, rho_ref)
136 class(unesco_eos), intent(in) :: this !< This EOS
137 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
138 real, intent(in) :: s !< Salinity [PSU]
139 real, intent(in) :: pressure !< Pressure [Pa]
140 real, intent(in) :: rho_ref !< A reference density [kg m-3]
141
142 ! Local variables
143 real :: t1 ! A copy of the temperature at a point [degC]
144 real :: s1 ! A copy of the salinity at a point [PSU]
145 real :: p1 ! Pressure converted to bars [bar]
146 real :: s12 ! The square root of salinity [PSU1/2]
147 real :: rho0 ! Density at 1 bar pressure [kg m-3]
148 real :: sig0 ! The anomaly of rho0 from R00 [kg m-3]
149 real :: ks ! The secant bulk modulus [bar]
150
151 p1 = pressure*1.0e-5 ; t1 = t
152 s1 = max(s, 0.0) ; s12 = sqrt(s1)
153
154 ! Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ).
155 sig0 = ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
156 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
157 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
158 rho0 = r00 + sig0
159
160 ! Compute rho(s,theta,p), first calculating the secant bulk modulus.
161 ks = (s000 + ( t1*(s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
162 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620))) )) + &
163 p1*( (s001 + ( t1*(s011 + t1*(s021 + t1*s031)) + &
164 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )) + &
165 p1*(s002 + ( t1*(s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )) )
166
167 density_anomaly_elem_unesco = ((r00 - rho_ref)*ks + (sig0*ks + p1*rho_ref)) / (ks - p1)
168
170
171!> In situ specific volume as fit by Jackett and McDougall, 1995 [m3 kg-1]
172!!
173!! This is an elemental function that can be applied to any combination of scalar and array inputs.
174real elemental function spec_vol_elem_unesco(this, t, s, pressure)
175 class(unesco_eos), intent(in) :: this !< This EOS
176 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
177 real, intent(in) :: s !< Salinity [PSU]
178 real, intent(in) :: pressure !< Pressure [Pa]
179
180 ! Local variables
181 real :: t1 ! A copy of the temperature at a point [degC]
182 real :: s1 ! A copy of the salinity at a point [PSU]
183 real :: p1 ! Pressure converted to bars [bar]
184 real :: s12 ! The square root of salinity [PSU1/2]l553
185 real :: rho0 ! Density at 1 bar pressure [kg m-3]
186 real :: ks ! The secant bulk modulus [bar]
187
188 p1 = pressure*1.0e-5 ; t1 = t
189 s1 = max(s, 0.0) ; s12 = sqrt(s1)
190
191 ! Compute rho(s,theta,p=0), which is the same as rho(s,t_insitu,p=0).
192 rho0 = r00 + ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
193 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
194 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
195
196 ! Compute rho(s,theta,p), first calculating the secant bulk modulus.
197 ks = (s000 + ( t1*(s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
198 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620))) )) + &
199 p1*( (s001 + ( t1*(s011 + t1*(s021 + t1*s031)) + &
200 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )) + &
201 p1*(s002 + ( t1*(s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )) )
202
203 spec_vol_elem_unesco = (ks - p1) / (rho0*ks)
204
205end function spec_vol_elem_unesco
206
207!> In situ specific volume anomaly as fit by Jackett and McDougall, 1995 [m3 kg-1]
208!!
209!! This is an elemental function that can be applied to any combination of scalar and array inputs.
210real elemental function spec_vol_anomaly_elem_unesco(this, t, s, pressure, spv_ref)
211 class(unesco_eos), intent(in) :: this !< This EOS
212 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
213 real, intent(in) :: s !< Salinity [PSU]
214 real, intent(in) :: pressure !< Pressure [Pa]
215 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
216
217 ! Local variables
218 real :: t1 ! A copy of the temperature at a point [degC]
219 real :: s1 ! A copy of the salinity at a point [PSU]
220 real :: p1 ! Pressure converted to bars [bar]
221 real :: s12 ! The square root of salinity [PSU1/2]
222 real :: rho0 ! Density at 1 bar pressure [kg m-3]
223 real :: ks ! The secant bulk modulus [bar]
224
225 p1 = pressure*1.0e-5 ; t1 = t
226 s1 = max(s, 0.0) ; s12 = sqrt(s1)
227
228 ! Compute rho(s,theta,p=0), which is the same as rho(s,t_insitu,p=0).
229 rho0 = r00 + ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
230 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
231 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
232
233 ! Compute rho(s,theta,p), first calculating the secant bulk modulus.
234 ks = (s000 + ( t1*(s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
235 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620))) )) + &
236 p1*( (s001 + ( t1*(s011 + t1*(s021 + t1*s031)) + &
237 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )) + &
238 p1*(s002 + ( t1*(s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )) )
239
240 spec_vol_anomaly_elem_unesco = (ks*(1.0 - (rho0*spv_ref)) - p1) / (rho0*ks)
241
243
244!> Calculate the partial derivatives of density with potential temperature and salinity
245!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995).
246elemental subroutine calculate_density_derivs_elem_unesco(this, T, S, pressure, drho_dT, drho_dS)
247 class(unesco_eos), intent(in) :: this !< This EOS
248 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
249 real, intent(in) :: s !< Salinity [PSU]
250 real, intent(in) :: pressure !< Pressure [Pa]
251 real, intent(out) :: drho_dt !< The partial derivative of density with potential
252 !! temperature [kg m-3 degC-1]
253 real, intent(out) :: drho_ds !< The partial derivative of density with salinity,
254 !! in [kg m-3 PSU-1]
255 ! Local variables
256 real :: t1 ! A copy of the temperature at a point [degC]
257 real :: s1 ! A copy of the salinity at a point [PSU]
258 real :: p1 ! Pressure converted to bars [bar]
259 real :: s12 ! The square root of salinity [PSU1/2]
260 real :: rho0 ! Density at 1 bar pressure [kg m-3]
261 real :: ks ! The secant bulk modulus [bar]
262 real :: drho0_dt ! Derivative of rho0 with T [kg m-3 degC-1]
263 real :: drho0_ds ! Derivative of rho0 with S [kg m-3 PSU-1]
264 real :: dks_dt ! Derivative of ks with T [bar degC-1]
265 real :: dks_ds ! Derivative of ks with S [bar psu-1]
266 real :: i_denom ! 1.0 / (ks - p1) [bar-1]
267
268 p1 = pressure*1.0e-5 ; t1 = t
269 s1 = max(s, 0.0) ; s12 = sqrt(s1)
270
271 ! Compute rho(s,theta,p=0) and its derivatives with temperature and salinity
272 rho0 = r00 + ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
273 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
274 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
275 drho0_dt = r01 + ( t1*(2.0*r02 + t1*(3.0*r03 + t1*(4.0*r04 + t1*(5.0*r05)))) + &
276 s1*(r11 + (t1*(2.0*r12 + t1*(3.0*r13 + t1*(4.0*r14))) + &
277 s12*(r61 + t1*(2.0*r62)) )) )
278 drho0_ds = r10 + ( t1*(r11 + t1*(r12 + t1*(r13 + t1*r14))) + &
279 (1.5*(s12*(r60 + t1*(r61 + t1*r62))) + s1*(2.0*r20)) )
280
281 ! Compute the secant bulk modulus and its derivatives with temperature and salinity
282 ks = ( s000 + (t1*(s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
283 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620)))) ) + &
284 p1*( (s001 + ( t1*(s011 + t1*(s021 + t1*s031)) + &
285 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )) + &
286 p1*(s002 + ( t1*(s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )) )
287 dks_dt = ( s010 + (t1*(2.0*s020 + t1*(3.0*s030 + t1*(4.0*s040))) + &
288 s1*((s110 + t1*(2.0*s120 + t1*(3.0*s130))) + s12*(s610 + t1*(2.0*s620)))) ) + &
289 p1*(((s011 + t1*(2.0*s021 + t1*(3.0*s031))) + s1*(s111 + t1*(2.0*s121)) ) + &
290 p1*(s012 + t1*(2.0*s022) + s1*(s112 + t1*(2.0*s122))) )
291 dks_ds = ( s100 + (t1*(s110 + t1*(s120 + t1*s130)) + 1.5*(s12*(s600 + t1*(s610 + t1*s620)))) ) + &
292 p1*((s101 + t1*(s111 + t1*s121) + s12*(1.5*s601)) + &
293 p1*(s102 + t1*(s112 + t1*s122)) )
294
295 i_denom = 1.0 / (ks - p1)
296 drho_dt = (ks*drho0_dt - dks_dt*((rho0*p1)*i_denom)) * i_denom
297 drho_ds = (ks*drho0_ds - dks_ds*((rho0*p1)*i_denom)) * i_denom
298
300
301!> Calculate second derivatives of density with respect to temperature, salinity, and pressure,
302!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995)
303elemental subroutine calculate_density_second_derivs_elem_unesco(this, T, S, pressure, &
304 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
305 class(unesco_eos), intent(in) :: this !< This EOS
306 real, intent(in) :: t !< Potential temperature referenced to 0 dbar [degC]
307 real, intent(in) :: s !< Salinity [PSU]
308 real, intent(in) :: pressure !< Pressure [Pa]
309 real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
310 !! to S [kg m-3 PSU-2]
311 real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
312 !! to T [kg m-3 PSU-1 degC-1]
313 real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
314 !! to T [kg m-3 degC-2]
315 real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
316 !! to pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1]
317 real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
318 !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
319
320 ! Local variables
321 real :: t1 ! A copy of the temperature at a point [degC]
322 real :: s1 ! A copy of the salinity at a point [PSU]
323 real :: p1 ! Pressure converted to bars [bar]
324 real :: s12 ! The square root of salinity [PSU1/2]
325 real :: i_s12 ! The inverse of the square root of salinity [PSU-1/2]
326 real :: rho0 ! Density at 1 bar pressure [kg m-3]
327 real :: drho0_dt ! Derivative of rho0 with T [kg m-3 degC-1]
328 real :: drho0_ds ! Derivative of rho0 with S [kg m-3 PSU-1]
329 real :: d2rho0_ds2 ! Second derivative of rho0 with salinity [kg m-3 PSU-1]
330 real :: d2rho0_dsdt ! Second derivative of rho0 with temperature and salinity [kg m-3 degC-1 PSU-1]
331 real :: d2rho0_dt2 ! Second derivative of rho0 with temperature [kg m-3 degC-2]
332 real :: ks ! The secant bulk modulus [bar]
333 real :: ks_0 ! The secant bulk modulus at zero pressure [bar]
334 real :: ks_1 ! The linear pressure dependence of the secant bulk modulus at zero pressure [nondim]
335 real :: ks_2 ! The quadratic pressure dependence of the secant bulk modulus at zero pressure [bar-1]
336 real :: dks_dp ! The derivative of the secant bulk modulus with pressure [nondim]
337 real :: dks_dt ! Derivative of the secant bulk modulus with temperature [bar degC-1]
338 real :: dks_ds ! Derivative of the secant bulk modulus with salinity [bar psu-1]
339 real :: d2ks_dt2 ! Second derivative of the secant bulk modulus with temperature [bar degC-2]
340 real :: d2ks_dsdt ! Second derivative of the secant bulk modulus with salinity and temperature [bar psu-1 degC-1]
341 real :: d2ks_ds2 ! Second derivative of the secant bulk modulus with salinity [bar psu-2]
342 real :: d2ks_dsdp ! Second derivative of the secant bulk modulus with salinity and pressure [psu-1]
343 real :: d2ks_dtdp ! Second derivative of the secant bulk modulus with temperature and pressure [degC-1]
344 real :: i_denom ! The inverse of the denominator of the expression for density [bar-1]
345
346 p1 = pressure*1.0e-5 ; t1 = t
347 s1 = max(s, 0.0) ; s12 = sqrt(s1)
348 ! The UNESCO equation of state is a fit to density, but it chooses a form that exhibits a
349 ! singularity in the second derivatives with salinity for fresh water. To avoid this, the
350 ! square root of salinity can be treated with a floor such that the contribution from the
351 ! S**1.5 terms to both the surface density and the secant bulk modulus are lost to roundoff.
352 ! This salinity is given by (~1e-16*S000/S600)**(2/3) ~= 3e-8 PSU, or S12 ~= 1.7e-4
353 i_s12 = 1.0 / (max(s12, 1.0e-4))
354
355 ! Calculate the density at sea level pressure and its derivatives
356 rho0 = r00 + ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
357 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
358 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
359 drho0_dt = r01 + ( t1*(2.0*r02 + t1*(3.0*r03 + t1*(4.0*r04 + t1*(5.0*r05)))) + &
360 s1*(r11 + ( t1*(2.0*r12 + t1*(3.0*r13 + t1*(4.0*r14))) + &
361 s12*(r61 + t1*(2.0*r62)) ) ) )
362 drho0_ds = r10 + ( t1*(r11 + t1*(r12 + t1*(r13 + t1*r14))) + &
363 (1.5*(s12*(r60 + t1*(r61 + t1*r62))) + s1*(2.0*r20)) )
364 d2rho0_ds2 = 0.75*(r60 + t1*(r61 + t1*r62))*i_s12 + 2.0*r20
365 d2rho0_dsdt = r11 + ( t1*(2.0*r12 + t1*(3.0*r13 + t1*(4.0*r14))) + s12*(1.5*r61 + t1*(3.0*r62)) )
366 d2rho0_dt2 = 2.0*r02 + ( t1*(6.0*r03 + t1*(12.0*r04 + t1*(20.0*r05))) + &
367 s1*((2.0*r12 + t1*(6.0*r13 + t1*(12.0*r14))) + s12*(2.0*r62)) )
368
369 ! Calculate the secant bulk modulus and its derivatives
370 ks_0 = s000 + ( t1*( s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
371 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620))) )
372 ks_1 = s001 + ( t1*( s011 + t1*(s021 + t1*s031)) + &
373 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )
374 ks_2 = s002 + ( t1*( s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )
375
376 ks = ks_0 + p1*(ks_1 + p1*ks_2)
377 dks_dp = ks_1 + 2.0*p1*ks_2
378 dks_dt = (s010 + ( t1*(2.0*s020 + t1*(3.0*s030 + t1*(4.0*s040))) + &
379 s1*((s110 + t1*(2.0*s120 + t1*(3.0*s130))) + s12*(s610 + t1*(2.0*s620))) )) + &
380 p1*((s011 + t1*(2.0*s021 + t1*(3.0*s031)) + s1*(s111 + t1*(2.0*s121))) + &
381 p1*(s012 + t1*(2.0*s022) + s1*(s112 + t1*(2.0*s122))))
382 dks_ds = (s100 + ( t1*(s110 + t1*(s120 + t1*s130)) + 1.5*(s12*(s600 + t1*(s610 + t1*s620))) )) + &
383 p1*((s101 + t1*(s111 + t1*s121) + s12*(1.5*s601)) + &
384 p1*(s102 + t1*(s112 + t1*s122)))
385 d2ks_ds2 = 0.75*((s600 + t1*(s610 + t1*s620)) + p1*s601)*i_s12
386 d2ks_dsdt = (s110 + ( t1*(2.0*s120 + t1*(3.0*s130)) + s12*(1.5*s610 + t1*(3.0*s620)) )) + &
387 p1*((s111 + t1*(2.0*s121)) + p1*(s112 + t1*(2.0*s122)))
388 d2ks_dt2 = 2.0*(s020 + ( t1*(3.0*s030 + t1*(6.0*s040)) + s1*((s120 + t1*(3.0*s130)) + s12*s620) )) + &
389 2.0*p1*((s021 + (t1*(3.0*s031) + s1*s121)) + p1*(s022 + s1*s122))
390
391 d2ks_dsdp = (s101 + (t1*(s111 + t1*s121) + s12*(1.5*s601))) + &
392 2.0*p1*(s102 + t1*(s112 + t1*s122))
393 d2ks_dtdp = (s011 + (t1*(2.0*s021 + t1*(3.0*s031)) + s1*(s111 + t1*(2.0*s121)))) + &
394 2.0*p1*(s012 + t1*(2.0*s022) + s1*(s112 + t1*(2.0*s122)))
395 i_denom = 1.0 / (ks - p1)
396
397 ! Expressions for density and its first derivatives are copied here for reference:
398 ! rho = rho0*ks * I_denom
399 ! drho_dT = I_denom*(ks*drho0_dT - p1*rho0*I_denom*dks_dT)
400 ! drho_dS = I_denom*(ks*drho0_dS - p1*rho0*I_denom*dks_dS)
401 ! drho_dp = 1.0e-5 * (rho0 * I_denom**2) * (ks - dks_dp*p1)
402
403 ! Finally calculate the second derivatives
404 drho_ds_ds = i_denom * ( ks*d2rho0_ds2 - (p1*i_denom) * &
405 (2.0*drho0_ds*dks_ds + rho0*(d2ks_ds2 - 2.0*dks_ds**2*i_denom)) )
406 drho_ds_dt = i_denom * (ks * d2rho0_dsdt - (p1*i_denom) * &
407 ((drho0_dt*dks_ds + drho0_ds*dks_dt) + &
408 rho0*(d2ks_dsdt - 2.0*(dks_ds*dks_dt)*i_denom)) )
409 drho_dt_dt = i_denom * ( ks*d2rho0_dt2 - (p1*i_denom) * &
410 (2.0*drho0_dt*dks_dt + rho0*(d2ks_dt2 - 2.0*dks_dt**2*i_denom)) )
411
412 ! The factor of 1.0e-5 is because pressure here is in bars, not Pa.
413 drho_ds_dp = (1.0e-5 * i_denom**2) * ( (ks*drho0_ds - rho0*dks_ds) - &
414 p1*( (dks_dp*drho0_ds + rho0*d2ks_dsdp) - &
415 2.0*(rho0*dks_ds) * ((dks_dp - 1.0)*i_denom) ) )
416 drho_dt_dp = (1.0e-5 * i_denom**2) * ( (ks*drho0_dt - rho0*dks_dt) - &
417 p1*( (dks_dp*drho0_dt + rho0*d2ks_dtdp) - &
418 2.0*(rho0*dks_dt) * ((dks_dp - 1.0)*i_denom) ) )
419
421
422!> Calculate the partial derivatives of specific volume with temperature and salinity
423!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995).
424elemental subroutine calculate_specvol_derivs_elem_unesco(this, T, S, pressure, dSV_dT, dSV_dS)
425 class(unesco_eos), intent(in) :: this !< This EOS
426 real, intent(in) :: t !< Potential temperature [degC]
427 real, intent(in) :: s !< Salinity [PSU]
428 real, intent(in) :: pressure !< Pressure [Pa]
429 real, intent(inout) :: dsv_dt !< The partial derivative of specific volume with
430 !! potential temperature [m3 kg-1 degC-1]
431 real, intent(inout) :: dsv_ds !< The partial derivative of specific volume with
432 !! salinity [m3 kg-1 PSU-1]
433 ! Local variables
434 real :: t1 ! A copy of the temperature at a point [degC]
435 real :: s1 ! A copy of the salinity at a point [PSU]
436 real :: p1 ! Pressure converted to bars [bar]
437 real :: s12 ! The square root of salinity [PSU1/2]
438 real :: rho0 ! Density at 1 bar pressure [kg m-3]
439 real :: ks ! The secant bulk modulus [bar]
440 real :: drho0_dt ! Derivative of rho0 with T [kg m-3 degC-1]
441 real :: drho0_ds ! Derivative of rho0 with S [kg m-3 PSU-1]
442 real :: dks_dt ! Derivative of ks with T [bar degC-1]
443 real :: dks_ds ! Derivative of ks with S [bar psu-1]
444 real :: i_denom2 ! 1.0 / (rho0*ks)**2 [m6 kg-2 bar-2]
445
446 p1 = pressure*1.0e-5 ; t1 = t
447 s1 = max(s, 0.0) ; s12 = sqrt(s1)
448
449 ! Compute rho(s,theta,p=0) and its derivatives with temperature and salinity
450 rho0 = r00 + ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
451 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
452 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
453 drho0_dt = r01 + ( t1*(2.0*r02 + t1*(3.0*r03 + t1*(4.0*r04 + t1*(5.0*r05)))) + &
454 s1*(r11 + (t1*(2.0*r12 + t1*(3.0*r13 + t1*(4.0*r14))) + &
455 s12*(r61 + t1*(2.0*r62)) )) )
456 drho0_ds = r10 + ( t1*(r11 + t1*(r12 + t1*(r13 + t1*r14))) + &
457 (1.5*(s12*(r60 + t1*(r61 + t1*r62))) + s1*(2.0*r20)) )
458
459 ! Compute the secant bulk modulus and its derivatives with temperature and salinity
460 ks = ( s000 + (t1*(s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
461 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620)))) ) + &
462 p1*( (s001 + ( t1*(s011 + t1*(s021 + t1*s031)) + &
463 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )) + &
464 p1*(s002 + ( t1*(s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )) )
465 dks_dt = ( s010 + (t1*(2.0*s020 + t1*(3.0*s030 + t1*(4.0*s040))) + &
466 s1*((s110 + t1*(2.0*s120 + t1*(3.0*s130))) + s12*(s610 + t1*(2.0*s620)))) ) + &
467 p1*(((s011 + t1*(2.0*s021 + t1*(3.0*s031))) + s1*(s111 + t1*(2.0*s121)) ) + &
468 p1*(s012 + t1*(2.0*s022) + s1*(s112 + t1*(2.0*s122))) )
469 dks_ds = ( s100 + (t1*(s110 + t1*(s120 + t1*s130)) + 1.5*(s12*(s600 + t1*(s610 + t1*s620)))) ) + &
470 p1*((s101 + t1*(s111 + t1*s121) + s12*(1.5*s601)) + &
471 p1*(s102 + t1*(s112 + t1*s122)) )
472
473 ! specvol = (ks - p1) / (rho0*ks) = 1/rho0 - p1/(rho0*ks)
474 i_denom2 = 1.0 / (rho0*ks)**2
475 dsv_dt = ((p1*rho0)*dks_dt + ((p1 - ks)*ks)*drho0_dt) * i_denom2
476 dsv_ds = ((p1*rho0)*dks_ds + ((p1 - ks)*ks)*drho0_ds) * i_denom2
477
479
480!> Compute the in situ density of sea water (rho) and the compressibility (drho/dp == C_sound^-2)
481!! at the given salinity, potential temperature and pressure using the UNESCO (1981)
482!! equation of state, as refit by Jackett and McDougall (1995).
483elemental subroutine calculate_compress_elem_unesco(this, T, S, pressure, rho, drho_dp)
484 class(unesco_eos), intent(in) :: this !< This EOS
485 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
486 real, intent(in) :: s !< Salinity [PSU]
487 real, intent(in) :: pressure !< Pressure [Pa]
488 real, intent(out) :: rho !< In situ density [kg m-3]
489 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
490 !! (also the inverse of the square of sound speed)
491 !! [s2 m-2]
492 ! Local variables
493 real :: t1 ! A copy of the temperature at a point [degC]
494 real :: s1 ! A copy of the salinity at a point [PSU]
495 real :: p1 ! Pressure converted to bars [bar]
496 real :: s12 ! The square root of salinity [PSU1/2]
497 real :: rho0 ! Density at 1 bar pressure [kg m-3]
498 real :: ks ! The secant bulk modulus [bar]
499 real :: ks_0 ! The secant bulk modulus at zero pressure [bar]
500 real :: ks_1 ! The linear pressure dependence of the secant bulk modulus at zero pressure [nondim]
501 real :: ks_2 ! The quadratic pressure dependence of the secant bulk modulus at zero pressure [bar-1]
502 real :: dks_dp ! The derivative of the secant bulk modulus with pressure [nondim]
503 real :: i_denom ! 1.0 / (ks - p1) [bar-1]
504
505 p1 = pressure*1.0e-5 ; t1 = t
506 s1 = max(s, 0.0) ; s12 = sqrt(s1)
507
508 ! Compute rho(s,theta,p=0), which is the same as rho(s,t_insitu,p=0).
509
510 rho0 = r00 + ( t1*(r01 + t1*(r02 + t1*(r03 + t1*(r04 + t1*r05)))) + &
511 s1*((r10 + t1*(r11 + t1*(r12 + t1*(r13 + t1*r14)))) + &
512 (s12*(r60 + t1*(r61 + t1*r62)) + s1*r20)) )
513
514 ! Calculate the secant bulk modulus and its derivative with pressure.
515 ks_0 = s000 + ( t1*( s010 + t1*(s020 + t1*(s030 + t1*s040))) + &
516 s1*((s100 + t1*(s110 + t1*(s120 + t1*s130))) + s12*(s600 + t1*(s610 + t1*s620))) )
517 ks_1 = s001 + ( t1*( s011 + t1*(s021 + t1*s031)) + &
518 s1*((s101 + t1*(s111 + t1*s121)) + s12*s601) )
519 ks_2 = s002 + ( t1*( s012 + t1*s022) + s1*(s102 + t1*(s112 + t1*s122)) )
520
521 ks = ks_0 + p1*(ks_1 + p1*ks_2)
522 dks_dp = ks_1 + 2.0*p1*ks_2
523 i_denom = 1.0 / (ks - p1)
524
525 ! Compute the in situ density, rho(s,theta,p), and its derivative with pressure.
526 rho = rho0*ks * i_denom
527 ! The factor of 1.0e-5 is because pressure here is in bars, not Pa.
528 drho_dp = 1.0e-5 * ((rho0 * (ks - p1*dks_dp)) * i_denom**2)
529
531
532!> Return the range of temperatures, salinities and pressures for which Jackett and McDougall (1995)
533!! refit the UNESCO equation of state has been fitted to observations. Care should be taken when
534!! applying this equation of state outside of its fit range.
535subroutine eos_fit_range_unesco(this, T_min, T_max, S_min, S_max, p_min, p_max)
536 class(unesco_eos), intent(in) :: this !< This EOS
537 real, optional, intent(out) :: T_min !< The minimum potential temperature over which this EoS is fitted [degC]
538 real, optional, intent(out) :: T_max !< The maximum potential temperature over which this EoS is fitted [degC]
539 real, optional, intent(out) :: S_min !< The minimum practical salinity over which this EoS is fitted [PSU]
540 real, optional, intent(out) :: S_max !< The maximum practical salinity over which this EoS is fitted [PSU]
541 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
542 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
543
544 if (present(t_min)) t_min = -2.5
545 if (present(t_max)) t_max = 40.0
546 if (present(s_min)) s_min = 0.0
547 if (present(s_max)) s_max = 42.0
548 if (present(p_min)) p_min = 0.0
549 if (present(p_max)) p_max = 1.0e8
550
551end subroutine eos_fit_range_unesco
552
553!> \namespace mom_eos_UNESCO
554!!
555!! \section section_EOS_UNESCO UNESCO (Jackett & McDougall) equation of state
556!!
557!! The UNESCO (1981) equation of state is an internationally defined standard fit valid over the
558!! range of pressures up to 10000 dbar, temperatures between the freezing point and 40 degC, and
559!! salinities between 0 and 42 PSU. Unfortunately, these expressions used in situ temperatures,
560!! whereas ocean models (including MOM6) effectively use potential temperatures as their state
561!! variables. To avoid needing multiple conversions, Jackett and McDougall (1995) refit the
562!! UNESCO equation of state to take potential temperature as a state variable, over the same
563!! valid range and functional form as the original UNESCO expressions. It is this refit from
564!! Jackett and McDougall (1995) that is coded up in this module.
565!!
566!! The functional form of the equation of state includes terms proportional to salinity to the
567!! 3/2 power. This introduces a singularity in the second derivative of density with salinity
568!! at a salinity of 0, but this has been addressed here by setting a floor of 1e-8 PSU on the
569!! salinity that is used in the denominator of these second derivative expressions. This value
570!! was chosen to imply a contribution that is smaller than numerical roundoff in the expression
571!! for density, which is the field for which the UNESCO equation of state was originally derived.
572!!
573!! Originally coded in 1999 by J. Stephens, revised in 2023 to unambiguously specify the order
574!! of arithmetic with parenthesis in every real sum of three or more terms.
575!!
576!! \subsection section_EOS_UNESCO_references References
577!!
578!! Gill, A. E., 1982: Atmosphere-Ocean Dynamics. Academic Press, 662 pp.
579!!
580!! Jackett, D. and T. McDougall, 1995: Minimal adjustment of hydrographic profiles to
581!! achieve static stability. J. Atmos. Ocean. Tech., 12, 381-389.
582!!
583!! UNESCO, 1981: Tenth report of the joint panel on oceanographic tables and standards.
584!! UNESCO Technical Papers in Marine Sci. No. 36, UNESCO, Paris.
585
586end module mom_eos_unesco