MOM_EOS_Jackett06.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 et al 2006 expressions that are often used in Hycom
7
9
10implicit none ; private
11
12public jackett06_eos
13
14!>@{ Parameters in the Jackett et al. equation of state, which is a fit to the Fiestel (2003)
15! equation of state for the range: -2 < theta < 40 [degC], 0 < S < 42 [PSU], 0 < p < 1e8 [Pa].
16! The notation here is for terms in the numerator of the expression for density of
17! RNabc for terms proportional to S**a * T**b * P**c, and terms in the denominator as RDabc.
18! For terms proportional to S**1.5, 6 is used in this notation.
19
20! --- coefficients for 25-term rational function sigloc().
21real, parameter :: &
22 rn000 = 9.9984085444849347d+02, & ! Density numerator constant coefficient [kg m-3]
23 rn001 = 1.1798263740430364d-06, & ! Density numerator P coefficient [kg m-3 Pa-1]
24 rn002 = -2.5862187075154352d-16, & ! Density numerator P^2 coefficient [kg m-3 Pa-2]
25 rn010 = 7.3471625860981584d+00, & ! Density numerator T coefficient [kg m-3 degC-1]
26 rn020 = -5.3211231792841769d-02, & ! Density numerator T^2 coefficient [kg m-3 degC-2]
27 rn021 = 9.8920219266399117d-12, & ! Density numerator T^2 P coefficient [kg m-3 degC-2 Pa-1]
28 rn022 = -3.2921414007960662d-20, & ! Density numerator T^2 P^2 coefficient [kg m-3 degC-2 Pa-2]
29 rn030 = 3.6492439109814549d-04, & ! Density numerator T^3 coefficient [kg m-3 degC-3]
30 rn100 = 2.5880571023991390d+00, & ! Density numerator S coefficient [kg m-3 PSU-1]
31 rn101 = 4.6996642771754730d-10, & ! Density numerator S P coefficient [kg m-3 PSU-1 Pa-1]
32 rn110 = -6.7168282786692355d-03, & ! Density numerator S T coefficient [kg m-3 degC-1 PSU-1]
33 rn200 = 1.9203202055760151d-03, & ! Density numerator S^2 coefficient [kg m-3]
34
35 rd001 = 6.7103246285651894d-10, & ! Density denominator P coefficient [Pa-1]
36 rd010 = 7.2815210113327091d-03, & ! Density denominator T coefficient [degC-1]
37 rd013 = -9.1534417604289062d-30, & ! Density denominator T P^3 coefficient [degC-1 Pa-3]
38 rd020 = -4.4787265461983921d-05, & ! Density denominator T^2 coefficient [degC-2]
39 rd030 = 3.3851002965802430d-07, & ! Density denominator T^3 coefficient [degC-3]
40 rd032 = -2.4461698007024582d-25, & ! Density denominator T^3 P^2 coefficient [degC-3 Pa-2]
41 rd040 = 1.3651202389758572d-10, & ! Density denominator T^4 coefficient [degC-4]
42 rd100 = 1.7632126669040377d-03, & ! Density denominator S coefficient [PSU-1]
43 rd110 = -8.8066583251206474d-06, & ! Density denominator S T coefficient [degC-1 PSU-1]
44 rd130 = -1.8832689434804897d-10, & ! Density denominator S T^3 coefficient [degC-3 PSU-1]
45 rd600 = 5.7463776745432097d-06, & ! Density denominator S^1.5 coefficient [PSU-1.5]
46 rd620 = 1.4716275472242334d-09 ! Density denominator S^1.5 T^2 coefficient [degC-2 PSU-1.5]
47!>@}
48
49!> The EOS_base implementation of the Jackett et al, 2006, equation of state
50type, extends (eos_base) :: jackett06_eos
51
52contains
53 !> Implementation of the in-situ density as an elemental function [kg m-3]
54 procedure :: density_elem => density_elem_jackett06
55 !> Implementation of the in-situ density anomaly as an elemental function [kg m-3]
56 procedure :: density_anomaly_elem => density_anomaly_elem_jackett06
57 !> Implementation of the in-situ specific volume as an elemental function [m3 kg-1]
58 procedure :: spec_vol_elem => spec_vol_elem_jackett06
59 !> Implementation of the in-situ specific volume anomaly as an elemental function [m3 kg-1]
60 procedure :: spec_vol_anomaly_elem => spec_vol_anomaly_elem_jackett06
61 !> Implementation of the calculation of derivatives of density
62 procedure :: calculate_density_derivs_elem => calculate_density_derivs_elem_jackett06
63 !> Implementation of the calculation of second derivatives of density
64 procedure :: calculate_density_second_derivs_elem => calculate_density_second_derivs_elem_jackett06
65 !> Implementation of the calculation of derivatives of specific volume
66 procedure :: calculate_specvol_derivs_elem => calculate_specvol_derivs_elem_jackett06
67 !> Implementation of the calculation of compressibility
68 procedure :: calculate_compress_elem => calculate_compress_elem_jackett06
69 !> Implementation of the range query function
70 procedure :: eos_fit_range => eos_fit_range_jackett06
71
72end type jackett06_eos
73
74contains
75
76!> In situ density of sea water using Jackett et al., 2006 [kg m-3]
77!!
78!! This is an elemental function that can be applied to any combination of scalar and array inputs.
79real elemental function density_elem_jackett06(this, t, s, pressure)
80 class(jackett06_eos), intent(in) :: this !< This EOS
81 real, intent(in) :: t !< Potential temperature relative to the surface [degC].
82 real, intent(in) :: s !< Salinity [PSU].
83 real, intent(in) :: pressure !< Pressure [Pa].
84
85 ! Local variables
86 real :: num_stp ! State dependent part of the numerator of the rational expresion
87 ! for density [kg m-3]
88 real :: den ! Denominator of the rational expresion for density [nondim]
89 real :: den_stp ! State dependent part of the denominator of the rational expresion
90 ! for density [nondim]
91 real :: i_den ! The inverse of the denominator of the rational expresion for density [nondim]
92 real :: t2 ! Temperature squared [degC2]
93 real :: s1_2 ! Limited square root of salinity [PSU1/2]
94
95 s1_2 = sqrt(max(0.0,s))
96 t2 = t*t
97
98 num_stp = (t*(rn010 + t*(rn020 + t*rn030)) + &
99 s*(rn100 + (t*rn110 + s*rn200)) ) + &
100 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022)))
101 den = 1.0 + ((t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
102 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
103 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013)) )
104 i_den = 1.0 / den
105
106 density_elem_jackett06 = (rn000 + num_stp)*i_den
107
108end function density_elem_jackett06
109
110!> In situ density anomaly of sea water using Jackett et al., 2006 [kg m-3]
111!!
112!! This is an elemental function that can be applied to any combination of scalar and array inputs.
113real elemental function density_anomaly_elem_jackett06(this, t, s, pressure, rho_ref)
114 class(jackett06_eos), intent(in) :: this !< This EOS
115 real, intent(in) :: t !< Potential temperature relative to the surface [degC].
116 real, intent(in) :: s !< Salinity [PSU].
117 real, intent(in) :: pressure !< Pressure [Pa].
118 real, intent(in) :: rho_ref !< A reference density [kg m-3].
119
120 ! Local variables
121 real :: num_stp ! State dependent part of the numerator of the rational expresion
122 ! for density [kg m-3]
123 real :: den ! Denominator of the rational expresion for density [nondim]
124 real :: den_stp ! State dependent part of the denominator of the rational expresion
125 ! for density [nondim]
126 real :: i_den ! The inverse of the denominator of the rational expresion for density [nondim]
127 real :: t2 ! Temperature squared [degC2]
128 real :: s1_2 ! Limited square root of salinity [PSU1/2]
129 real :: rho0 ! The surface density of fresh water at 0 degC, perhaps less the refernce density [kg m-3]
130
131 s1_2 = sqrt(max(0.0,s))
132 t2 = t*t
133
134 num_stp = (t*(rn010 + t*(rn020 + t*rn030)) + &
135 s*(rn100 + (t*rn110 + s*rn200)) ) + &
136 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022)))
137 den = 1.0 + ((t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
138 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
139 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013)) )
140 i_den = 1.0 / den
141
142 rho0 = rn000 - rho_ref*den
143
144 density_anomaly_elem_jackett06 = (rho0 + num_stp)*i_den
145
147
148!> In situ specific volume of sea water using Jackett et al., 2006 [m3 kg-1]
149!!
150!! This is an elemental function that can be applied to any combination of scalar and array inputs.
151real elemental function spec_vol_elem_jackett06(this, t, s, pressure)
152 class(jackett06_eos), intent(in) :: this !< This EOS
153 real, intent(in) :: t !< potential temperature relative to the surface [degC].
154 real, intent(in) :: s !< salinity [PSU].
155 real, intent(in) :: pressure !< pressure [Pa].
156
157 ! Local variables
158 real :: num_stp ! State dependent part of the numerator of the rational expresion
159 ! for density (not specific volume) [kg m-3]
160 real :: den_stp ! State dependent part of the denominator of the rational expresion
161 ! for density (not specific volume) [nondim]
162 real :: i_num ! The inverse of the numerator of the rational expresion for density [nondim]
163 real :: t2 ! Temperature squared [degC2]
164 real :: s1_2 ! Limited square root of salinity [PSU1/2]
165
166 s1_2 = sqrt(max(0.0,s))
167 t2 = t*t
168
169 num_stp = (t*(rn010 + t*(rn020 + t*rn030)) + &
170 s*(rn100 + (t*rn110 + s*rn200)) ) + &
171 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022)))
172 den_stp = (t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
173 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
174 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013))
175 i_num = 1.0 / (rn000 + num_stp)
176
177 spec_vol_elem_jackett06 = (1.0 + den_stp) * i_num
178
179end function spec_vol_elem_jackett06
180
181!> In situ specific volume anomaly of sea water using Jackett et al., 2006 [m3 kg-1]
182!!
183!! This is an elemental function that can be applied to any combination of scalar and array inputs.
184real elemental function spec_vol_anomaly_elem_jackett06(this, t, s, pressure, spv_ref)
185 class(jackett06_eos), intent(in) :: this !< This EOS
186 real, intent(in) :: t !< potential temperature relative to the surface [degC].
187 real, intent(in) :: s !< salinity [PSU].
188 real, intent(in) :: pressure !< pressure [Pa].
189 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
190
191 ! Local variables
192 real :: num_stp ! State dependent part of the numerator of the rational expresion
193 ! for density (not specific volume) [kg m-3]
194 real :: den_stp ! State dependent part of the denominator of the rational expresion
195 ! for density (not specific volume) [nondim]
196 real :: i_num ! The inverse of the numerator of the rational expresion for density [nondim]
197 real :: t2 ! Temperature squared [degC2]
198 real :: s1_2 ! Limited square root of salinity [PSU1/2]
199
200 s1_2 = sqrt(max(0.0,s))
201 t2 = t*t
202
203 num_stp = (t*(rn010 + t*(rn020 + t*rn030)) + &
204 s*(rn100 + (t*rn110 + s*rn200)) ) + &
205 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022)))
206 den_stp = (t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
207 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
208 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013))
209 i_num = 1.0 / (rn000 + num_stp)
210
211 ! This form is slightly more complicated, but it cancels the leading terms better.
212 spec_vol_anomaly_elem_jackett06 = ((1.0 - spv_ref*rn000) + (den_stp - spv_ref*num_stp)) * i_num
213
215
216!> Calculate the partial derivatives of density with potential temperature and salinity
217!! using Jackett et al., 2006
218elemental subroutine calculate_density_derivs_elem_jackett06(this, T, S, pressure, drho_dT, drho_dS)
219 class(jackett06_eos), intent(in) :: this !< This EOS
220 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
221 real, intent(in) :: s !< Salinity [PSU]
222 real, intent(in) :: pressure !< Pressure [Pa]
223 real, intent(out) :: drho_dt !< The partial derivative of density with potential
224 !! temperature [kg m-3 degC-1]
225 real, intent(out) :: drho_ds !< The partial derivative of density with salinity,
226 !! in [kg m-3 PSU-1]
227 ! Local variables
228 real :: num ! Numerator of the rational expresion for density [kg m-3]
229 real :: den ! Denominator of the rational expresion for density [nondim]
230 real :: i_denom2 ! The inverse of the square of the denominator of the rational expression
231 ! for density [nondim]
232 real :: dnum_dt ! The derivative of num with potential temperature [kg m-3 degC-1]
233 real :: dnum_ds ! The derivative of num with salinity [kg m-3 PSU-1]
234 real :: dden_dt ! The derivative of den with potential temperature [degC-1]
235 real :: dden_ds ! The derivative of den with salinity PSU-1]
236 real :: t2 ! Temperature squared [degC2]
237 real :: s1_2 ! Limited square root of salinity [PSU1/2]
238
239 s1_2 = sqrt(max(0.0,s))
240 t2 = t*t
241
242 num = rn000 + ((t*(rn010 + t*(rn020 + t*rn030)) + &
243 s*(rn100 + (t*rn110 + s*rn200)) ) + &
244 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022))) )
245 den = 1.0 + ((t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
246 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
247 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013)) )
248
249 dnum_dt = ((rn010 + t*(2.*rn020 + t*(3.*rn030))) + s*rn110) + &
250 pressure*t*(2.*rn021 + pressure*(2.*rn022))
251 dnum_ds = (rn100 + (t*rn110 + s*(2.*rn200))) + pressure*rn101
252 dden_dt = ((rd010 + t*((2.*rd020) + t*((3.*rd030) + t*(4.*rd040)))) + &
253 s*((rd110 + t2*(3.*rd130)) + s1_2*t*(2.*rd620)) ) + &
254 pressure**2*(t2*3.*rd032 + pressure*rd013)
255 dden_ds = rd100 + (t*(rd110 + t2*rd130) + s1_2*(1.5*rd600 + t2*(1.5*rd620)))
256 i_denom2 = 1.0 / den**2
257
258 ! rho = num / den
259 drho_dt = (dnum_dt * den - num * dden_dt) * i_denom2
260 drho_ds = (dnum_ds * den - num * dden_ds) * i_denom2
261
263
264!> Calculate second derivatives of density with respect to temperature, salinity, and pressure,
265!! using Jackett et al., 2006
266elemental subroutine calculate_density_second_derivs_elem_jackett06(this, T, S, pressure, &
267 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
268 class(jackett06_eos), intent(in) :: this !< This EOS
269 real, intent(in) :: t !< Potential temperature referenced to 0 dbar [degC]
270 real, intent(in) :: s !< Salinity [PSU]
271 real, intent(in) :: pressure !< Pressure [Pa]
272 real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
273 !! to S [kg m-3 PSU-2]
274 real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
275 !! to T [kg m-3 PSU-1 degC-1]
276 real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
277 !! to T [kg m-3 degC-2]
278 real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
279 !! to pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1]
280 real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
281 !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
282
283 ! Local variables
284 real :: num ! Numerator of the rational expresion for density [kg m-3]
285 real :: den ! Denominator of the rational expresion for density [nondim]
286 real :: i_num2 ! The inverse of the square of the numerator of the rational expression
287 ! for density [nondim]
288 real :: dnum_dt ! The derivative of num with potential temperature [kg m-3 degC-1]
289 real :: dnum_ds ! The derivative of num with salinity [kg m-3 PSU-1]
290 real :: dden_dt ! The derivative of den with potential temperature [degC-1]
291 real :: dden_ds ! The derivative of den with salinity PSU-1]
292 real :: dnum_dp ! The derivative of num with pressure [kg m-3 dbar-1]
293 real :: dden_dp ! The derivative of det with pressure [dbar-1]
294 real :: d2num_dt2 ! The second derivative of num with potential temperature [kg m-3 degC-2]
295 real :: d2num_dt_ds ! The second derivative of num with potential temperature and
296 ! salinity [kg m-3 degC-1 PSU-1]
297 real :: d2num_ds2 ! The second derivative of num with salinity [kg m-3 PSU-2]
298 real :: d2num_dt_dp ! The second derivative of num with potential temperature and
299 ! pressure [kg m-3 degC-1 dbar-1]
300 real :: d2num_ds_dp ! The second derivative of num with salinity and
301 ! pressure [kg m-3 PSU-1 dbar-1]
302 real :: d2den_dt2 ! The second derivative of den with potential temperature [degC-2]
303 real :: d2den_dt_ds ! The second derivative of den with potential temperature and salinity [degC-1 PSU-1]
304 real :: d2den_ds2 ! The second derivative of den with salinity [PSU-2]
305 real :: d2den_dt_dp ! The second derivative of den with potential temperature and pressure [degC-1 dbar-1]
306 real :: d2den_ds_dp ! The second derivative of den with salinity and pressure [PSU-1 dbar-1]
307 real :: t2 ! Temperature squared [degC2]
308 real :: s1_2 ! Limited square root of salinity [PSU1/2]
309 real :: i_s12 ! The inverse of the square root of salinity [PSU-1/2]
310 real :: i_denom2 ! The inverse of the square of the denominator of the rational expression
311 ! for density [nondim]
312 real :: i_denom3 ! The inverse of the cube of the denominator of the rational expression
313 ! for density [nondim]
314
315 s1_2 = sqrt(max(0.0,s))
316 t2 = t*t
317
318 num = rn000 + ((t*(rn010 + t*(rn020 + t*rn030)) + &
319 s*(rn100 + (t*rn110 + s*rn200)) ) + &
320 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022))) )
321 den = 1.0 + ((t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
322 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
323 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013)) )
324 ! rho = num*I_den
325
326 dnum_dt = ((rn010 + t*(2.*rn020 + t*(3.*rn030))) + s*rn110) + &
327 pressure*t*(2.*rn021 + pressure*(2.*rn022))
328 dnum_ds = (rn100 + (t*rn110 + s*(2.*rn200))) + pressure*rn101
329 dnum_dp = rn001 + ((t2*rn021 + s*rn101) + pressure*(2.*rn002 + t2*(2.*rn022)))
330 d2num_dt2 = 2.*rn020 + t*(6.*rn030) + pressure*(2.*rn021 + pressure*(2.*rn022))
331 d2num_dt_ds = rn110
332 d2num_ds2 = 2.*rn200
333 d2num_dt_dp = t*(2.*rn021 + pressure*(4.*rn022))
334 d2num_ds_dp = rn101
335
336 dden_dt = ((rd010 + t*((2.*rd020) + t*((3.*rd030) + t*(4.*rd040)))) + &
337 s*((rd110 + t2*(3.*rd130)) + s1_2*t*(2.*rd620)) ) + &
338 pressure**2*(t2*3.*rd032 + pressure*rd013)
339 dden_ds = rd100 + (t*(rd110 + t2*rd130) + s1_2*(1.5*rd600 + t2*(1.5*rd620)))
340 dden_dp = rd001 + pressure*t*(t2*(2.*rd032) + pressure*(3.*rd013))
341
342 d2den_dt2 = (((2.*rd020) + t*((6.*rd030) + t*(12.*rd040))) + &
343 s*(t*(6.*rd130) + s1_2*(2.*rd620)) ) + pressure**2*(t*(6.*rd032))
344 d2den_dt_ds = (rd110 + t2*3.*rd130) + (t*s1_2)*(3.0*rd620)
345 d2den_dt_dp = pressure*(t2*(6.*rd032) + pressure*(3.*rd013))
346 d2den_ds_dp = 0.0
347
348 ! The Jackett et al. 2006 equation of state is a fit to density, but it chooses a form that
349 ! exhibits a singularity in the second derivatives with salinity for fresh water. To avoid
350 ! this, the square root of salinity can be treated with a floor such that the contribution from
351 ! the 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/RD600)**(2/3) ~= 7e-8 PSU, or S1_2 ~= 2.6e-4
353 i_s12 = 1.0 / (max(s1_2, 1.0e-4))
354 d2den_ds2 = (0.75*rd600 + t2*(0.75*rd620)) * i_s12
355
356 i_denom3 = 1.0 / den**3
357
358 ! In deriving the following, it is useful to note that:
359 ! drho_dp = (dnum_dp * den - num * dden_dp) / den**2
360 ! drho_dT = (dnum_dT * den - num * dden_dT) / den**2
361 ! drho_dS = (dnum_dS * den - num * dden_dS) / den**2
362 drho_ds_ds = (den*(den*d2num_ds2 - 2.*dnum_ds*dden_ds) + num*(2.*dden_ds**2 - den*d2den_ds2)) * i_denom3
363 drho_ds_dt = (den*(den*d2num_dt_ds - (dnum_dt*dden_ds + dnum_ds*dden_dt)) + &
364 num*(2.*dden_dt*dden_ds - den*d2den_dt_ds)) * i_denom3
365 drho_dt_dt = (den*(den*d2num_dt2 - 2.*dnum_dt*dden_dt) + num*(2.*dden_dt**2 - den*d2den_dt2)) * i_denom3
366
367 drho_ds_dp = (den*(den*d2num_ds_dp - (dnum_dp*dden_ds + dnum_ds*dden_dp)) + &
368 num*(2.*dden_ds*dden_dp - den*d2den_ds_dp)) * i_denom3
369 drho_dt_dp = (den*(den*d2num_dt_dp - (dnum_dp*dden_dt + dnum_dt*dden_dp)) + &
370 num*(2.*dden_dt*dden_dp - den*d2den_dt_dp)) * i_denom3
371
373
374!> Calculate second derivatives of density with respect to temperature, salinity, and pressure,
375!! using Jackett et al., 2006
376elemental subroutine calculate_specvol_derivs_elem_jackett06(this, T, S, pressure, dSV_dT, dSV_dS)
377 class(jackett06_eos), intent(in) :: this !< This EOS
378 real, intent(in) :: t !< Potential temperature [degC]
379 real, intent(in) :: s !< Salinity [PSU]
380 real, intent(in) :: pressure !< Pressure [Pa]
381 real, intent(inout) :: dsv_dt !< The partial derivative of specific volume with
382 !! potential temperature [m3 kg-1 degC-1]
383 real, intent(inout) :: dsv_ds !< The partial derivative of specific volume with
384 !! salinity [m3 kg-1 PSU-1]
385
386 ! Local variables
387 real :: num ! Numerator of the rational expresion for density (not specific volume) [kg m-3]
388 real :: den ! Denominator of the rational expresion for density (not specific volume) [nondim]
389 real :: i_num2 ! The inverse of the square of the numerator of the rational expression
390 ! for density [nondim]
391 real :: dnum_dt ! The derivative of num with potential temperature [kg m-3 degC-1]
392 real :: dnum_ds ! The derivative of num with salinity [kg m-3 PSU-1]
393 real :: dden_dt ! The derivative of den with potential temperature [degC-1]
394 real :: dden_ds ! The derivative of den with salinity PSU-1]
395 real :: t2 ! Temperature squared [degC2]
396 real :: s1_2 ! Limited square root of salinity [PSU1/2]
397
398 s1_2 = sqrt(max(0.0,s))
399 t2 = t*t
400
401 num = rn000 + ((t*(rn010 + t*(rn020 + t*rn030)) + &
402 s*(rn100 + (t*rn110 + s*rn200)) ) + &
403 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022))) )
404 den = 1.0 + ((t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
405 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
406 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013)) )
407
408 dnum_dt = ((rn010 + t*(2.*rn020 + t*(3.*rn030))) + s*rn110) + &
409 pressure*t*(2.*rn021 + pressure*(2.*rn022))
410 dnum_ds = (rn100 + (t*rn110 + s*(2.*rn200))) + pressure*rn101
411 dden_dt = ((rd010 + t*((2.*rd020) + t*((3.*rd030) + t*(4.*rd040)))) + &
412 s*((rd110 + t2*(3.*rd130)) + s1_2*t*(2.*rd620)) ) + &
413 pressure**2*(t2*3.*rd032 + pressure*rd013)
414 dden_ds = rd100 + (t*(rd110 + t2*rd130) + s1_2*(1.5*rd600 + t2*(1.5*rd620)))
415 i_num2 = 1.0 / num**2
416
417 ! SV = den / num
418 dsv_dt = (num * dden_dt - dnum_dt * den) * i_num2
419 dsv_ds = (num * dden_ds - dnum_ds * den) * i_num2
420
422
423!> Compute the in situ density of sea water (rho) and the compressibility (drho/dp == C_sound^-2)
424!! at the given salinity, potential temperature and pressure using Jackett et al., 2006
425elemental subroutine calculate_compress_elem_jackett06(this, T, S, pressure, rho, drho_dp)
426 class(jackett06_eos), intent(in) :: this !< This EOS
427 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
428 real, intent(in) :: s !< Salinity [PSU]
429 real, intent(in) :: pressure !< Pressure [Pa]
430 real, intent(out) :: rho !< In situ density [kg m-3]
431 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
432 !! (also the inverse of the square of sound speed)
433 !! [s2 m-2]
434 ! Local variables
435 real :: num ! Numerator of the rational expresion for density [kg m-3]
436 real :: den ! Denominator of the rational expresion for density [nondim]
437 real :: i_den ! The inverse of the denominator of the rational expression for density [nondim]
438 real :: dnum_dp ! The derivative of num with pressure [kg m-3 dbar-1]
439 real :: dden_dp ! The derivative of den with pressure [dbar-1]
440 real :: t2 ! Temperature squared [degC2]
441 real :: s1_2 ! Limited square root of salinity [PSU1/2]
442 integer :: j
443
444 s1_2 = sqrt(max(0.0,s))
445 t2 = t*t
446
447 num = rn000 + ((t*(rn010 + t*(rn020 + t*rn030)) + &
448 s*(rn100 + (t*rn110 + s*rn200)) ) + &
449 pressure*(rn001 + ((t2*rn021 + s*rn101) + pressure*(rn002 + t2*rn022))) )
450 den = 1.0 + ((t*(rd010 + t*(rd020 + t*(rd030 + t* rd040))) + &
451 s*(rd100 + (t*(rd110 + t2*rd130) + s1_2*(rd600 + t2*rd620))) ) + &
452 pressure*(rd001 + pressure*t*(t2*rd032 + pressure*rd013)) )
453 dnum_dp = rn001 + ((t2*rn021 + s*rn101) + pressure*(2.*rn002 + t2*(2.*rn022)))
454 dden_dp = rd001 + pressure*t*(t2*(2.*rd032) + pressure*(3.*rd013))
455
456 i_den = 1.0 / den
457 rho = num * i_den
458 drho_dp = (dnum_dp * den - num * dden_dp) * i_den**2
459
461
462!> Return the range of temperatures, salinities and pressures for which the Jackett et al. (2006)
463!! equation of state has been fitted to observations. Care should be taken when applying this
464!! equation of state outside of its fit range.
465subroutine eos_fit_range_jackett06(this, T_min, T_max, S_min, S_max, p_min, p_max)
466 class(jackett06_eos), intent(in) :: this !< This EOS
467 real, optional, intent(out) :: T_min !< The minimum potential temperature over which this EoS is fitted [degC]
468 real, optional, intent(out) :: T_max !< The maximum potential temperature over which this EoS is fitted [degC]
469 real, optional, intent(out) :: S_min !< The minimum practical salinity over which this EoS is fitted [PSU]
470 real, optional, intent(out) :: S_max !< The maximum practical salinity over which this EoS is fitted [PSU]
471 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
472 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
473
474 ! Note that the actual fit range is given for the surface range of temperatures and salinities,
475 ! but Jackett et al. use a more limited range of properties at higher pressures.
476 if (present(t_min)) t_min = -4.5
477 if (present(t_max)) t_max = 40.0
478 if (present(s_min)) s_min = 0.0
479 if (present(s_max)) s_max = 42.0
480 if (present(p_min)) p_min = 0.0
481 if (present(p_max)) p_max = 8.5e7
482
483end subroutine eos_fit_range_jackett06
484
485
486!> \namespace mom_eos_Jackett06
487!!
488!! \section section_EOS_Jackett06 Jackett et al. 2006 (Hycom-25-term) equation of state
489!!
490!! Jackett et al. (2006) provide an approximation for the in situ density as a function of
491!! potential temperature, salinity, and pressure. This 25 term equation of state is
492!! frequently used in Hycom for a potential density, at which point it only has 17 terms
493!! and so is commonly called the "17-term equation of state" there. Here the full expressions
494!! for the in situ densities are used.
495!!
496!! The functional form of this equation of state includes terms proportional to salinity to the
497!! 3/2 power. This introduces a singularity in the second derivative of density with salinity
498!! at a salinity of 0, but this has been addressed here by setting a floor of 1e-8 PSU on the
499!! salinity that is used in the denominator of these second derivative expressions. This value
500!! was chosen to imply a contribution that is smaller than numerical roundoff in the expression for
501!! density, which is the field for which the Jackett et al. equation of state was originally derived.
502!!
503!! \subsection section_EOS_Jackett06_references References
504!!
505!! Jackett, D., T. McDougall, R. Feistel, D. Wright and S. Griffies (2006),
506!! Algorithms for density, potential temperature, conservative
507!! temperature, and the freezing temperature of seawater, JAOT
508!! doi.org/10.1175/JTECH1946.1
509
510end module mom_eos_jackett06