MOM_EOS_Roquet_rho.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 expressions of Roquet et al. (2015) that are used in NEMO
7
9
10implicit none ; private
11
12public roquet_rho_eos
13
14real, parameter :: pa2kb = 1.e-8 !< Conversion factor between Pa and kbar [kbar Pa-1]
15!>@{ Parameters in the Roquet_rho (Roquet density) equation of state
16real, parameter :: rdeltas = 32. ! An offset to salinity before taking its square root [g kg-1]
17real, parameter :: r1_s0 = 0.875/35.16504 ! The inverse of a plausible range of oceanic salinities [kg g-1]
18real, parameter :: i_ts = 0.025 ! The inverse of a plausible range of oceanic temperatures [degC-1]
19
20! The following are the coefficients of the fit to the reference density profile (rho00p) as a function of
21! pressure (P), with a contribution R0c * P**(c+1). The nomenclature follows Roquet.
22real, parameter :: r00 = 4.6494977072e+01*pa2kb ! rho00p P coef. [kg m-3 Pa-1]
23real, parameter :: r01 = -5.2099962525*pa2kb**2 ! rho00p P**2 coef. [kg m-3 Pa-2]
24real, parameter :: r02 = 2.2601900708e-01*pa2kb**3 ! rho00p P**3 coef. [kg m-3 Pa-3]
25real, parameter :: r03 = 6.4326772569e-02*pa2kb**4 ! rho00p P**4 coef. [kg m-3 Pa-4]
26real, parameter :: r04 = 1.5616995503e-02*pa2kb**5 ! rho00p P**5 coef. [kg m-3 Pa-5]
27real, parameter :: r05 = -1.7243708991e-03*pa2kb**6 ! rho00p P**6 coef. [kg m-3 Pa-6]
28
29! The following are coefficients of contributions to density as a function of the square root
30! of normalized salinity with an offset (zs), temperature (T) and pressure (P), with a contribution
31! EOSabc * zs**a * T**b * P**c. The numbers here are copied directly from Roquet et al. (2015), but
32! the expressions here do not use the same nondimensionalization for pressure or temperature as they do.
33real, parameter :: eos000 = 8.0189615746e+02 ! A constant density contribution [kg m-3]
34real, parameter :: eos100 = 8.6672408165e+02 ! EoS zs coef. [kg m-3]
35real, parameter :: eos200 = -1.7864682637e+03 ! EoS zs**2 coef. [kg m-3]
36real, parameter :: eos300 = 2.0375295546e+03 ! EoS zs**3 coef. [kg m-3]
37real, parameter :: eos400 = -1.2849161071e+03 ! EoS zs**4 coef. [kg m-3]
38real, parameter :: eos500 = 4.3227585684e+02 ! EoS zs**5 coef. [kg m-3]
39real, parameter :: eos600 = -6.0579916612e+01 ! EoS zs**6 coef. [kg m-3]
40real, parameter :: eos010 = 2.6010145068e+01*i_ts ! EoS T coef. [kg m-3 degC-1]
41real, parameter :: eos110 = -6.5281885265e+01*i_ts ! EoS zs * T coef. [kg m-3 degC-1]
42real, parameter :: eos210 = 8.1770425108e+01*i_ts ! EoS zs**2 * T coef. [kg m-3 degC-1]
43real, parameter :: eos310 = -5.6888046321e+01*i_ts ! EoS zs**3 * T coef. [kg m-3 degC-1]
44real, parameter :: eos410 = 1.7681814114e+01*i_ts ! EoS zs**2 * T coef. [kg m-3 degC-1]
45real, parameter :: eos510 = -1.9193502195*i_ts ! EoS zs**5 * T coef. [kg m-3 degC-1]
46real, parameter :: eos020 = -3.7074170417e+01*i_ts**2 ! EoS T**2 coef. [kg m-3 degC-2]
47real, parameter :: eos120 = 6.1548258127e+01*i_ts**2 ! EoS zs * T**2 coef. [kg m-3 degC-2]
48real, parameter :: eos220 = -6.0362551501e+01*i_ts**2 ! EoS zs**2 * T**2 coef. [kg m-3 degC-2]
49real, parameter :: eos320 = 2.9130021253e+01*i_ts**2 ! EoS zs**3 * T**2 coef. [kg m-3 degC-2]
50real, parameter :: eos420 = -5.4723692739*i_ts**2 ! EoS zs**4 * T**2 coef. [kg m-3 degC-2]
51real, parameter :: eos030 = 2.1661789529e+01*i_ts**3 ! EoS T**3 coef. [kg m-3 degC-3]
52real, parameter :: eos130 = -3.3449108469e+01*i_ts**3 ! EoS zs * T**3 coef. [kg m-3 degC-3]
53real, parameter :: eos230 = 1.9717078466e+01*i_ts**3 ! EoS zs**2 * T**3 coef. [kg m-3 degC-3]
54real, parameter :: eos330 = -3.1742946532*i_ts**3 ! EoS zs**3 * T**3 coef. [kg m-3 degC-3]
55real, parameter :: eos040 = -8.3627885467*i_ts**4 ! EoS T**4 coef. [kg m-3 degC-4]
56real, parameter :: eos140 = 1.1311538584e+01*i_ts**4 ! EoS zs * T**4 coef. [kg m-3 degC-4]
57real, parameter :: eos240 = -5.3563304045*i_ts**4 ! EoS zs**2 * T**4 coef. [kg m-3 degC-4]
58real, parameter :: eos050 = 5.4048723791e-01*i_ts**5 ! EoS T**5 coef. [kg m-3 degC-5]
59real, parameter :: eos150 = 4.8169980163e-01*i_ts**5 ! EoS zs * T**5 coef. [kg m-3 degC-5]
60real, parameter :: eos060 = -1.9083568888e-01*i_ts**6 ! EoS T**6 [kg m-3 degC-6]
61real, parameter :: eos001 = 1.9681925209e+01*pa2kb ! EoS P coef. [kg m-3 Pa-1]
62real, parameter :: eos101 = -4.2549998214e+01*pa2kb ! EoS zs * P coef. [kg m-3 Pa-1]
63real, parameter :: eos201 = 5.0774768218e+01*pa2kb ! EoS zs**2 * P coef. [kg m-3 Pa-1]
64real, parameter :: eos301 = -3.0938076334e+01*pa2kb ! EoS zs**3 * P coef. [kg m-3 Pa-1]
65real, parameter :: eos401 = 6.6051753097*pa2kb ! EoS zs**4 * P coef. [kg m-3 Pa-1]
66real, parameter :: eos011 = -1.3336301113e+01*(i_ts*pa2kb) ! EoS T * P coef. [kg m-3 degC-1 Pa-1]
67real, parameter :: eos111 = -4.4870114575*(i_ts*pa2kb) ! EoS zs * T * P coef. [kg m-3 degC-1 Pa-1]
68real, parameter :: eos211 = 5.0042598061*(i_ts*pa2kb) ! EoS zs**2 * T * P coef. [kg m-3 degC-1 Pa-1]
69real, parameter :: eos311 = -6.5399043664e-01*(i_ts*pa2kb) ! EoS zs**3 * T * P coef. [kg m-3 degC-1 Pa-1]
70real, parameter :: eos021 = 6.7080479603*(i_ts**2*pa2kb) ! EoS T**2 * P coef. [kg m-3 degC-2 Pa-1]
71real, parameter :: eos121 = 3.5063081279*(i_ts**2*pa2kb) ! EoS zs * T**2 * P coef. [kg m-3 degC-2 Pa-1]
72real, parameter :: eos221 = -1.8795372996*(i_ts**2*pa2kb) ! EoS zs**2 * T**2 * P coef. [kg m-3 degC-2 Pa-1]
73real, parameter :: eos031 = -2.4649669534*(i_ts**3*pa2kb) ! EoS T**3 * P coef. [kg m-3 degC-3 Pa-1]
74real, parameter :: eos131 = -5.5077101279e-01*(i_ts**3*pa2kb) ! EoS zs * T**3 * P coef. [kg m-3 degC-3 Pa-1]
75real, parameter :: eos041 = 5.5927935970e-01*(i_ts**4*pa2kb) ! EoS T**4 * P coef. [kg m-3 degC-4 Pa-1]
76real, parameter :: eos002 = 2.0660924175*pa2kb**2 ! EoS P**2 coef. [kg m-3 Pa-2]
77real, parameter :: eos102 = -4.9527603989*pa2kb**2 ! EoS zs * P**2 coef. [kg m-3 Pa-2]
78real, parameter :: eos202 = 2.5019633244*pa2kb**2 ! EoS zs**2 * P**2 coef. [kg m-3 Pa-2]
79real, parameter :: eos012 = 2.0564311499*(i_ts*pa2kb**2) ! EoS T * P**2 coef. [kg m-3 degC-1 Pa-2]
80real, parameter :: eos112 = -2.1311365518e-01*(i_ts*pa2kb**2) ! EoS zs * T * P**2 coef. [kg m-3 degC-1 Pa-2]
81real, parameter :: eos022 = -1.2419983026*(i_ts**2*pa2kb**2) ! EoS T**2 * P**2 coef. [kg m-3 degC-2 Pa-2]
82real, parameter :: eos003 = -2.3342758797e-02*pa2kb**3 ! EoS P**3 coef. [kg m-3 Pa-3]
83real, parameter :: eos103 = -1.8507636718e-02*pa2kb**3 ! EoS zs * P**3 coef. [kg m-3 Pa-3]
84real, parameter :: eos013 = 3.7969820455e-01*(i_ts*pa2kb**3) ! EoS T * P**3 coef. [kg m-3 degC-1 Pa-3]
85
86real, parameter :: alp000 = eos010 ! Constant in the drho_dT fit [kg m-3 degC-1]
87real, parameter :: alp100 = eos110 ! drho_dT fit zs coef. [kg m-3 degC-1]
88real, parameter :: alp200 = eos210 ! drho_dT fit zs**2 coef. [kg m-3 degC-1]
89real, parameter :: alp300 = eos310 ! drho_dT fit zs**3 coef. [kg m-3 degC-1]
90real, parameter :: alp400 = eos410 ! drho_dT fit zs**4 coef. [kg m-3 degC-1]
91real, parameter :: alp500 = eos510 ! drho_dT fit zs**5 coef. [kg m-3 degC-1]
92real, parameter :: alp010 = 2.*eos020 ! drho_dT fit T coef. [kg m-3 degC-2]
93real, parameter :: alp110 = 2.*eos120 ! drho_dT fit zs * T coef. [kg m-3 degC-2]
94real, parameter :: alp210 = 2.*eos220 ! drho_dT fit zs**2 * T coef. [kg m-3 degC-2]
95real, parameter :: alp310 = 2.*eos320 ! drho_dT fit zs**3 * T coef. [kg m-3 degC-2]
96real, parameter :: alp410 = 2.*eos420 ! drho_dT fit zs**4 * T coef. [kg m-3 degC-2]
97real, parameter :: alp020 = 3.*eos030 ! drho_dT fit T**2 coef. [kg m-3 degC-3]
98real, parameter :: alp120 = 3.*eos130 ! drho_dT fit zs * T**2 coef. [kg m-3 degC-3]
99real, parameter :: alp220 = 3.*eos230 ! drho_dT fit zs**2 * T**2 coef. [kg m-3 degC-3]
100real, parameter :: alp320 = 3.*eos330 ! drho_dT fit zs**3 * T**2 coef. [kg m-3 degC-3]
101real, parameter :: alp030 = 4.*eos040 ! drho_dT fit T**3 coef. [kg m-3 degC-4]
102real, parameter :: alp130 = 4.*eos140 ! drho_dT fit zs * T**3 coef. [kg m-3 degC-4]
103real, parameter :: alp230 = 4.*eos240 ! drho_dT fit zs**2 * T**3 coef. [kg m-3 degC-4]
104real, parameter :: alp040 = 5.*eos050 ! drho_dT fit T**4 coef. [kg m-3 degC-5]
105real, parameter :: alp140 = 5.*eos150 ! drho_dT fit zs* * T**4 coef. [kg m-3 degC-5]
106real, parameter :: alp050 = 6.*eos060 ! drho_dT fit T**5 coef. [kg m-3 degC-6]
107real, parameter :: alp001 = eos011 ! drho_dT fit P coef. [kg m-3 degC-1 Pa-1]
108real, parameter :: alp101 = eos111 ! drho_dT fit zs * P coef. [kg m-3 degC-1 Pa-1]
109real, parameter :: alp201 = eos211 ! drho_dT fit zs**2 * P coef. [kg m-3 degC-1 Pa-1]
110real, parameter :: alp301 = eos311 ! drho_dT fit zs**3 * P coef. [kg m-3 degC-1 Pa-1]
111real, parameter :: alp011 = 2.*eos021 ! drho_dT fit T * P coef. [kg m-3 degC-2 Pa-1]
112real, parameter :: alp111 = 2.*eos121 ! drho_dT fit zs * T * P coef. [kg m-3 degC-2 Pa-1]
113real, parameter :: alp211 = 2.*eos221 ! drho_dT fit zs**2 * T * P coef. [kg m-3 degC-2 Pa-1]
114real, parameter :: alp021 = 3.*eos031 ! drho_dT fit T**2 * P coef. [kg m-3 degC-3 Pa-1]
115real, parameter :: alp121 = 3.*eos131 ! drho_dT fit zs * T**2 * P coef. [kg m-3 degC-3 Pa-1]
116real, parameter :: alp031 = 4.*eos041 ! drho_dT fit T**3 * P coef. [kg m-3 degC-4 Pa-1]
117real, parameter :: alp002 = eos012 ! drho_dT fit P**2 coef. [kg m-3 degC-1 Pa-2]
118real, parameter :: alp102 = eos112 ! drho_dT fit zs * P**2 coef. [kg m-3 degC-1 Pa-2]
119real, parameter :: alp012 = 2.*eos022 ! drho_dT fit T * P**2 coef. [kg m-3 degC-2 Pa-2]
120real, parameter :: alp003 = eos013 ! drho_dT fit P**3 coef. [kg m-3 degC-1 Pa-3]
121
122real, parameter :: bet000 = 0.5*eos100*r1_s0 ! Constant in the drho_dS fit [kg m-3 ppt-1]
123real, parameter :: bet100 = eos200*r1_s0 ! drho_dS fit zs coef. [kg m-3 ppt-1]
124real, parameter :: bet200 = 1.5*eos300*r1_s0 ! drho_dS fit zs**2 coef. [kg m-3 ppt-1]
125real, parameter :: bet300 = 2.0*eos400*r1_s0 ! drho_dS fit zs**3 coef. [kg m-3 ppt-1]
126real, parameter :: bet400 = 2.5*eos500*r1_s0 ! drho_dS fit zs**4 coef. [kg m-3 ppt-1]
127real, parameter :: bet500 = 3.0*eos600*r1_s0 ! drho_dS fit zs**5 coef. [kg m-3 ppt-1]
128real, parameter :: bet010 = 0.5*eos110*r1_s0 ! drho_dS fit T coef. [kg m-3 ppt-1 degC-1]
129real, parameter :: bet110 = eos210*r1_s0 ! drho_dS fit zs * T coef. [kg m-3 ppt-1 degC-1]
130real, parameter :: bet210 = 1.5*eos310*r1_s0 ! drho_dS fit zs**2 * T coef. [kg m-3 ppt-1 degC-1]
131real, parameter :: bet310 = 2.0*eos410*r1_s0 ! drho_dS fit zs**3 * T coef. [kg m-3 ppt-1 degC-1]
132real, parameter :: bet410 = 2.5*eos510*r1_s0 ! drho_dS fit zs**4 * T coef. [kg m-3 ppt-1 degC-1]
133real, parameter :: bet020 = 0.5*eos120*r1_s0 ! drho_dS fit T**2 coef. [kg m-3 ppt-1 degC-2]
134real, parameter :: bet120 = eos220*r1_s0 ! drho_dS fit zs * T**2 coef. [kg m-3 ppt-1 degC-2]
135real, parameter :: bet220 = 1.5*eos320*r1_s0 ! drho_dS fit zs**2 * T**2 coef. [kg m-3 ppt-1 degC-2]
136real, parameter :: bet320 = 2.0*eos420*r1_s0 ! drho_dS fit zs**3 * T**2 coef. [kg m-3 ppt-1 degC-2]
137real, parameter :: bet030 = 0.5*eos130*r1_s0 ! drho_dS fit T**3 coef. [kg m-3 ppt-1 degC-3]
138real, parameter :: bet130 = eos230*r1_s0 ! drho_dS fit zs * T**3 coef. [kg m-3 ppt-1 degC-3]
139real, parameter :: bet230 = 1.5*eos330*r1_s0 ! drho_dS fit zs**2 * T**3 coef. [kg m-3 ppt-1 degC-3]
140real, parameter :: bet040 = 0.5*eos140*r1_s0 ! drho_dS fit T**4 coef. [kg m-3 ppt-1 degC-4]
141real, parameter :: bet140 = eos240*r1_s0 ! drho_dS fit zs * T**4 coef. [kg m-3 ppt-1 degC-4]
142real, parameter :: bet050 = 0.5*eos150*r1_s0 ! drho_dS fit T**5 coef. [kg m-3 ppt-1 degC-5]
143real, parameter :: bet001 = 0.5*eos101*r1_s0 ! drho_dS fit P coef. [kg m-3 ppt-1 Pa-1]
144real, parameter :: bet101 = eos201*r1_s0 ! drho_dS fit zs * P coef. [kg m-3 ppt-1 Pa-1]
145real, parameter :: bet201 = 1.5*eos301*r1_s0 ! drho_dS fit zs**2 * P coef. [kg m-3 ppt-1 Pa-1]
146real, parameter :: bet301 = 2.0*eos401*r1_s0 ! drho_dS fit zs**3 * P coef. [kg m-3 ppt-1 Pa-1]
147real, parameter :: bet011 = 0.5*eos111*r1_s0 ! drho_dS fit T * P coef. [kg m-3 ppt-1 degC-1 Pa-1]
148real, parameter :: bet111 = eos211*r1_s0 ! drho_dS fit zs * T * P coef. [kg m-3 ppt-1 degC-1 Pa-1]
149real, parameter :: bet211 = 1.5*eos311*r1_s0 ! drho_dS fit zs**2 * T * P coef. [kg m-3 ppt-1 degC-1 Pa-1]
150real, parameter :: bet021 = 0.5*eos121*r1_s0 ! drho_dS fit T**2 * P coef. [kg m-3 ppt-1 degC-2 Pa-1]
151real, parameter :: bet121 = eos221*r1_s0 ! drho_dS fit zs * T**2 * P coef. [kg m-3 ppt-1 degC-2 Pa-1]
152real, parameter :: bet031 = 0.5*eos131*r1_s0 ! drho_dS fit T**3 * P coef. [kg m-3 ppt-1 degC-3 Pa-1]
153real, parameter :: bet002 = 0.5*eos102*r1_s0 ! drho_dS fit P**2 coef. [kg m-3 ppt-1 Pa-2]
154real, parameter :: bet102 = eos202*r1_s0 ! drho_dS fit zs * P**2 coef. [kg m-3 ppt-1 Pa-2]
155real, parameter :: bet012 = 0.5*eos112*r1_s0 ! drho_dS fit T * P**2 coef. [kg m-3 ppt-1 degC-1 Pa-2]
156real, parameter :: bet003 = 0.5*eos103*r1_s0 ! drho_dS fit P**3 coef. [kg m-3 ppt-1 Pa-3]
157!>@}
158
159!> The EOS_base implementation of the Roquet et al., 2015, equation of state
160type, extends (eos_base) :: roquet_rho_eos
161
162contains
163 !> Implementation of the in-situ density as an elemental function [kg m-3]
164 procedure :: density_elem => density_elem_roquet_rho
165 !> Implementation of the in-situ density anomaly as an elemental function [kg m-3]
166 procedure :: density_anomaly_elem => density_anomaly_elem_roquet_rho
167 !> Implementation of the in-situ specific volume as an elemental function [m3 kg-1]
168 procedure :: spec_vol_elem => spec_vol_elem_roquet_rho
169 !> Implementation of the in-situ specific volume anomaly as an elemental function [m3 kg-1]
170 procedure :: spec_vol_anomaly_elem => spec_vol_anomaly_elem_roquet_rho
171 !> Implementation of the calculation of derivatives of density
172 procedure :: calculate_density_derivs_elem => calculate_density_derivs_elem_roquet_rho
173 !> Implementation of the calculation of second derivatives of density
174 procedure :: calculate_density_second_derivs_elem => calculate_density_second_derivs_elem_roquet_rho
175 !> Implementation of the calculation of derivatives of specific volume
176 procedure :: calculate_specvol_derivs_elem => calculate_specvol_derivs_elem_roquet_rho
177 !> Implementation of the calculation of compressibility
178 procedure :: calculate_compress_elem => calculate_compress_elem_roquet_rho
179 !> Implementation of the range query function
180 procedure :: eos_fit_range => eos_fit_range_roquet_rho
181
182 !> Local implementation of generic calculate_density_array for efficiency
183 procedure :: calculate_density_array => calculate_density_array_roquet_rho
184 !> Local implementation of generic calculate_spec_vol_array for efficiency
185 procedure :: calculate_spec_vol_array => calculate_spec_vol_array_roquet_rho
186
187end type roquet_rho_eos
188
189contains
190
191!> In situ density of sea water from Roquet et al., 2015 [kg m-3]
192!!
193!! This is an elemental function that can be applied to any combination of scalar and array inputs.
194real elemental function density_elem_roquet_rho(this, t, s, pressure)
195 class(roquet_rho_eos), intent(in) :: this !< This EOS
196 real, intent(in) :: t !< Conservative temperature [degC]
197 real, intent(in) :: s !< Absolute salinity [g kg-1]
198 real, intent(in) :: pressure !< Pressure [Pa]
199
200 ! Local variables
201 real :: zp ! Pressure [Pa]
202 real :: zt ! Conservative temperature [degC]
203 real :: zs ! The square root of absolute salinity with an offset normalized
204 ! by an assumed salinity range [nondim]
205 real :: rho00p ! A pressure-dependent but temperature and salinity independent contribution to
206 ! density at the reference temperature and salinity [kg m-3]
207 real :: rhots ! Density without a pressure-dependent contribution [kg m-3]
208 real :: rhots0 ! A contribution to density from temperature and salinity anomalies at the
209 ! surface pressure [kg m-3]
210 real :: rhots1 ! A density contribution proportional to pressure [kg m-3 Pa-1]
211 real :: rhots2 ! A density contribution proportional to pressure**2 [kg m-3 Pa-2]
212 real :: rhots3 ! A density contribution proportional to pressure**3 [kg m-3 Pa-3]
213 real :: rho0s0 ! Salinity dependent density at the surface pressure and zero temperature [kg m-3]
214
215 ! The following algorithm was published by Roquet et al. (2015), intended for use with NEMO.
216
217 ! Conversions to the units used here.
218 zt = t
219 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
220 zp = pressure
221
222 ! The next two lines should be used if it is necessary to convert potential temperature and
223 ! practical salinity to conservative temperature and absolute salinity.
224 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
225 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
226
227 rhots3 = eos003 + (zs*eos103 + zt*eos013)
228 rhots2 = eos002 + (zs*(eos102 + zs*eos202) &
229 + zt*(eos012 + (zs*eos112 + zt*eos022)) )
230 rhots1 = eos001 + (zs*(eos101 + zs*(eos201 + zs*(eos301 + zs*eos401))) &
231 + zt*(eos011 + (zs*(eos111 + zs*(eos211 + zs*eos311)) &
232 + zt*(eos021 + (zs*(eos121 + zs*eos221) &
233 + zt*(eos031 + (zs*eos131 + zt*eos041)) )) )) )
234 rhots0 = zt*(eos010 &
235 + (zs*(eos110 + zs*(eos210 + zs*(eos310 + zs*(eos410 + zs*eos510)))) &
236 + zt*(eos020 + (zs*(eos120 + zs*(eos220 + zs*(eos320 + zs*eos420))) &
237 + zt*(eos030 + (zs*(eos130 + zs*(eos230 + zs*eos330)) &
238 + zt*(eos040 + (zs*(eos140 + zs*eos240) &
239 + zt*(eos050 + (zs*eos150 + zt*eos060)) )) )) )) ) )
240
241 rho0s0 = eos000 + zs*(eos100 + zs*(eos200 + zs*(eos300 + zs*(eos400 + zs*(eos500 + zs*eos600)))))
242
243 rho00p = zp*(r00 + zp*(r01 + zp*(r02 + zp*(r03 + zp*(r04 + zp*r05)))))
244
245 rhots = (rhots0 + rho0s0) + zp*(rhots1 + zp*(rhots2 + zp*rhots3))
246 density_elem_roquet_rho = rhots + rho00p ! In situ density [kg m-3]
247
248end function density_elem_roquet_rho
249
250!> In situ density anomaly of sea water from Roquet et al., 2015 [kg m-3]
251!!
252!! This is an elemental function that can be applied to any combination of scalar and array inputs.
253real elemental function density_anomaly_elem_roquet_rho(this, t, s, pressure, rho_ref)
254 class(roquet_rho_eos), intent(in) :: this !< This EOS
255 real, intent(in) :: t !< Conservative temperature [degC]
256 real, intent(in) :: s !< Absolute salinity [g kg-1]
257 real, intent(in) :: pressure !< Pressure [Pa]
258 real, intent(in) :: rho_ref !< A reference density [kg m-3]
259
260 ! Local variables
261 real :: zp ! Pressure [Pa]
262 real :: zt ! Conservative temperature [degC]
263 real :: zs ! The square root of absolute salinity with an offset normalized
264 ! by an assumed salinity range [nondim]
265 real :: rho00p ! A pressure-dependent but temperature and salinity independent contribution to
266 ! density at the reference temperature and salinity [kg m-3]
267 real :: rhots ! Density without a pressure-dependent contribution [kg m-3]
268 real :: rhots0 ! A contribution to density from temperature and salinity anomalies at the
269 ! surface pressure [kg m-3]
270 real :: rhots1 ! A density contribution proportional to pressure [kg m-3 Pa-1]
271 real :: rhots2 ! A density contribution proportional to pressure**2 [kg m-3 Pa-2]
272 real :: rhots3 ! A density contribution proportional to pressure**3 [kg m-3 Pa-3]
273 real :: rho0s0 ! Salinity dependent density at the surface pressure and zero temperature [kg m-3]
274
275 ! The following algorithm was published by Roquet et al. (2015), intended for use with NEMO.
276
277 ! Conversions to the units used here.
278 zt = t
279 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
280 zp = pressure
281
282 ! The next two lines should be used if it is necessary to convert potential temperature and
283 ! practical salinity to conservative temperature and absolute salinity.
284 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
285 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
286
287 rhots3 = eos003 + (zs*eos103 + zt*eos013)
288 rhots2 = eos002 + (zs*(eos102 + zs*eos202) &
289 + zt*(eos012 + (zs*eos112 + zt*eos022)) )
290 rhots1 = eos001 + (zs*(eos101 + zs*(eos201 + zs*(eos301 + zs*eos401))) &
291 + zt*(eos011 + (zs*(eos111 + zs*(eos211 + zs*eos311)) &
292 + zt*(eos021 + (zs*(eos121 + zs*eos221) &
293 + zt*(eos031 + (zs*eos131 + zt*eos041)) )) )) )
294 rhots0 = zt*(eos010 &
295 + (zs*(eos110 + zs*(eos210 + zs*(eos310 + zs*(eos410 + zs*eos510)))) &
296 + zt*(eos020 + (zs*(eos120 + zs*(eos220 + zs*(eos320 + zs*eos420))) &
297 + zt*(eos030 + (zs*(eos130 + zs*(eos230 + zs*eos330)) &
298 + zt*(eos040 + (zs*(eos140 + zs*eos240) &
299 + zt*(eos050 + (zs*eos150 + zt*eos060)) )) )) )) ) )
300
301 rho0s0 = eos000 + zs*(eos100 + zs*(eos200 + zs*(eos300 + zs*(eos400 + zs*(eos500 + zs*eos600)))))
302
303 rho00p = zp*(r00 + zp*(r01 + zp*(r02 + zp*(r03 + zp*(r04 + zp*r05)))))
304
305 rho0s0 = rho0s0 - rho_ref
306
307 rhots = (rhots0 + rho0s0) + zp*(rhots1 + zp*(rhots2 + zp*rhots3))
308 density_anomaly_elem_roquet_rho = rhots + rho00p ! In situ density [kg m-3]
309
311
312!> In situ specific volume of sea water from Roquet et al., 2015 [kg m-3]
313!!
314!! This is an elemental function that can be applied to any combination of scalar and array inputs.
315real elemental function spec_vol_elem_roquet_rho(this, t, s, pressure)
316 class(roquet_rho_eos), intent(in) :: this !< This EOS
317 real, intent(in) :: t !< Conservative temperature [degC]
318 real, intent(in) :: s !< Absolute salinity [g kg-1]
319 real, intent(in) :: pressure !< Pressure [Pa]
320
321 spec_vol_elem_roquet_rho = 1. / density_elem_roquet_rho(this, t, s, pressure)
322
323end function spec_vol_elem_roquet_rho
324
325!> In situ specific volume anomaly of sea water from Roquet et al., 2015 [kg m-3]
326!!
327!! This is an elemental function that can be applied to any combination of scalar and array inputs.
328real elemental function spec_vol_anomaly_elem_roquet_rho(this, t, s, pressure, spv_ref)
329 class(roquet_rho_eos), intent(in) :: this !< This EOS
330 real, intent(in) :: t !< Conservative temperature [degC]
331 real, intent(in) :: s !< Absolute salinity [g kg-1]
332 real, intent(in) :: pressure !< Pressure [Pa]
333 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
334
337
339
340!> For a given thermodynamic state, calculate the derivatives of density with conservative
341!! temperature and absolute salinity, using the density polynomial fit EOS from Roquet et al. (2015).
342elemental subroutine calculate_density_derivs_elem_roquet_rho(this, T, S, pressure, drho_dT, drho_dS)
343 class(roquet_rho_eos), intent(in) :: this !< This EOS
344 real, intent(in) :: t !< Conservative temperature [degC]
345 real, intent(in) :: s !< Absolute salinity [g kg-1]
346 real, intent(in) :: pressure !< Pressure [Pa]
347 real, intent(out) :: drho_dt !< The partial derivative of density with potential
348 !! temperature [kg m-3 degC-1]
349 real, intent(out) :: drho_ds !< The partial derivative of density with salinity,
350 !! in [kg m-3 ppt-1]
351
352 ! Local variables
353 real :: zp ! Pressure [Pa]
354 real :: zt ! Conservative temperature [degC]
355 real :: zs ! The square root of absolute salinity with an offset normalized
356 ! by an assumed salinity range [nondim]
357 real :: drdzt0 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1]
358 ! from temperature anomalies at the surface pressure
359 real :: drdzt1 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1 Pa-1]
360 ! proportional to pressure
361 real :: drdzt2 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1 Pa-2]
362 ! proportional to pressure**2
363 real :: drdzt3 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1 Pa-3]
364 ! proportional to pressure**3
365 real :: drdzs0 ! A contribution to the partial derivative of density with
366 ! salinity [kg m-3 ppt-1] from temperature anomalies at the surface pressure
367 real :: drdzs1 ! A contribution to the partial derivative of density with
368 ! salinity [kg m-3 ppt-1 Pa-1] proportional to pressure
369 real :: drdzs2 ! A contribution to the partial derivative of density with
370 ! salinity [kg m-3 ppt-1 Pa-2] proportional to pressure**2
371 real :: drdzs3 ! A contribution to the partial derivative of density with
372 ! salinity [kg m-3 ppt-1 Pa-3] proportional to pressure**3
373
374 ! Conversions to the units used here.
375 zt = t
376 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
377 zp = pressure
378
379 ! The next two lines should be used if it is necessary to convert potential temperature and
380 ! practical salinity to conservative temperature and absolute salinity.
381 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
382 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
383
384 ! Find the partial derivative of density with temperature
385 drdzt3 = alp003
386 drdzt2 = alp002 + (zs*alp102 + zt*alp012)
387 drdzt1 = alp001 + (zs*(alp101 + zs*(alp201 + zs*alp301)) &
388 + zt*(alp011 + (zs*(alp111 + zs*alp211) &
389 + zt*(alp021 + (zs*alp121 + zt*alp031)) )) )
390 drdzt0 = alp000 + (zs*(alp100 + zs*(alp200 + zs*(alp300 + zs*(alp400 + zs*alp500)))) &
391 + zt*(alp010 + (zs*(alp110 + zs*(alp210 + zs*(alp310 + zs*alp410))) &
392 + zt*(alp020 + (zs*(alp120 + zs*(alp220 + zs*alp320)) &
393 + zt*(alp030 + (zt*(alp040 + (zs*alp140 + zt*alp050)) &
394 + zs*(alp130 + zs*alp230) )) )) )) )
395
396 drho_dt = drdzt0 + zp*(drdzt1 + zp*(drdzt2 + zp*drdzt3))
397
398 ! Find the partial derivative of density with salinity
399 drdzs3 = bet003
400 drdzs2 = bet002 + (zs*bet102 + zt*bet012)
401 drdzs1 = bet001 + (zs*(bet101 + zs*(bet201 + zs*bet301)) &
402 + zt*(bet011 + (zs*(bet111 + zs*bet211) &
403 + zt*(bet021 + (zs*bet121 + zt*bet031)) )) )
404 drdzs0 = bet000 + (zs*(bet100 + zs*(bet200 + zs*(bet300 + zs*(bet400 + zs*bet500)))) &
405 + zt*(bet010 + (zs*(bet110 + zs*(bet210 + zs*(bet310 + zs*bet410))) &
406 + zt*(bet020 + (zs*(bet120 + zs*(bet220 + zs*bet320)) &
407 + zt*(bet030 + (zt*(bet040 + (zs*bet140 + zt*bet050)) &
408 + zs*(bet130 + zs*bet230) )) )) )) )
409
410 ! The division by zs here is because zs = sqrt(S + S0), so drho_dS = dzs_dS * drho_dzs = (0.5 / zs) * drho_dzs
411 drho_ds = (drdzs0 + zp*(drdzs1 + zp*(drdzs2 + zp * drdzs3))) / zs
412
414
415!> Second derivatives of density with respect to temperature, salinity, and pressure
416elemental subroutine calculate_density_second_derivs_elem_roquet_rho(this, T, S, pressure, &
417 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
418 class(roquet_rho_eos), intent(in) :: this !< This EOS
419 real, intent(in) :: t !< Conservative temperature [degC]
420 real, intent(in) :: s !< Absolute salinity [g kg-1]
421 real, intent(in) :: pressure !< Pressure [Pa]
422 real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
423 !! to S [kg m-3 ppt-2]
424 real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
425 !! to T [kg m-3 ppt-1 degC-1]
426 real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
427 !! to T [kg m-3 degC-2]
428 real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
429 !! to pressure [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1]
430 real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
431 !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
432
433 ! Local variables
434 real :: zp ! Pressure [Pa]
435 real :: zt ! Conservative temperature [degC]
436 real :: zs ! The square root of absolute salinity with an offset normalized
437 ! by an assumed salinity range [nondim]
438 real :: i_s ! The inverse of zs [nondim]
439 real :: d2r_p0 ! A contribution to one of the second derivatives that is independent of pressure [various]
440 real :: d2r_p1 ! A contribution to one of the second derivatives that is proportional to pressure [various]
441 real :: d2r_p2 ! A contribution to one of the second derivatives that is proportional to pressure**2 [various]
442 real :: d2r_p3 ! A contribution to one of the second derivatives that is proportional to pressure**3 [various]
443
444 ! Conversions to the units used here.
445 zt = t
446 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
447 zp = pressure
448
449 ! The next two lines should be used if it is necessary to convert potential temperature and
450 ! practical salinity to conservative temperature and absolute salinity.
451 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
452 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
453
454 i_s = 1.0 / zs
455
456 ! Find drho_ds_ds
457 d2r_p3 = -eos103*i_s**2
458 d2r_p2 = -(eos102 + zt*eos112)*i_s**2
459 d2r_p1 = (3.*eos301 + (zt*(3.*eos311) + zs*(8.*eos401))) &
460 - ( eos101 + zt*(eos111 + zt*(eos121 + zt*eos131)) )*i_s**2
461 d2r_p0 = (3.*eos300 + (zs*(8.*eos400 + zs*(15.*eos500 + zs*(24.*eos600))) &
462 + zt*(3.*eos310 + (zs*(8.*eos410 + zs*(15.*eos510)) &
463 + zt*(3.*eos320 + (zs*(8.*eos420) + zt*(3.*eos330))) )) )) &
464 - (eos100 + zt*(eos110 + zt*(eos120 + zt*(eos130 + zt*(eos140 + zt*eos150)))) )*i_s**2
465 drho_ds_ds = (0.5*r1_s0)**2 * ((d2r_p0 + zp*(d2r_p1 + zp*(d2r_p2 + zp*d2r_p3))) * i_s)
466
467 ! Find drho_ds_dt
468 d2r_p2 = eos112
469 d2r_p1 = eos111 + (zs*(2.*eos211 + zs*(3.*eos311)) &
470 + zt*(2.*eos121 + (zs*(4.*eos221) + zt*(3.*eos131))) )
471 d2r_p0 = eos110 + (zs*(2.*eos210 + zs*(3.*eos310 + zs*(4.*eos410 + zs*(5.*eos510)))) &
472 + zt*(2.*eos120 + (zs*(4.*eos220 + zs*(6.*eos320 + zs*(8.*eos420))) &
473 + zt*(3.*eos130 + (zs*(6.*eos230 + zs*(9.*eos330)) &
474 + zt*(4.*eos140 + (zs*(8.*eos240) &
475 + zt*(5.*eos150))) )) )) )
476 drho_ds_dt = (0.5*r1_s0) * ((d2r_p0 + zp*(d2r_p1 + zp*d2r_p2)) * i_s)
477
478 ! Find drho_dt_dt
479 d2r_p2 = 2.*eos022
480 d2r_p1 = 2.*eos021 + (zs*(2.*eos121 + zs*(2.*eos221)) &
481 + zt*(6.*eos031 + (zs*(6.*eos131) + zt*(12.*eos041))) )
482 d2r_p0 = 2.*eos020 + (zs*(2.*eos120 + zs*( 2.*eos220 + zs*( 2.*eos320 + zs * (2.*eos420)))) &
483 + zt*(6.*eos030 + (zs*( 6.*eos130 + zs*( 6.*eos230 + zs * (6.*eos330))) &
484 + zt*(12.*eos040 + (zs*(12.*eos140 + zs *(12.*eos240)) &
485 + zt*(20.*eos050 + (zs*(20.*eos150) &
486 + zt*(30.*eos060) )) )) )) )
487 drho_dt_dt = (d2r_p0 + zp*(d2r_p1 + zp*d2r_p2))
488
489 ! Find drho_ds_dp
490 d2r_p2 = 3.*eos103
491 d2r_p1 = 2.*eos102 + (zs*(4.*eos202) + zt*(2.*eos112))
492 d2r_p0 = eos101 + (zs*(2.*eos201 + zs*(3.*eos301 + zs*(4.*eos401))) &
493 + zt*(eos111 + (zs*(2.*eos211 + zs*(3.*eos311)) &
494 + zt*( eos121 + (zs*(2.*eos221) + zt*eos131)) )) )
495 drho_ds_dp = ((d2r_p0 + zp*(d2r_p1 + zp*d2r_p2)) * i_s) * (0.5*r1_s0)
496
497 ! Find drho_dt_dp
498 d2r_p2 = 3.*eos013
499 d2r_p1 = 2.*eos012 + (zs*(2.*eos112) + zt*(4.*eos022))
500 d2r_p0 = eos011 + (zs*(eos111 + zs*( eos211 + zs* eos311)) &
501 + zt*(2.*eos021 + (zs*(2.*eos121 + zs*(2.*eos221)) &
502 + zt*(3.*eos031 + (zs*(3.*eos131) + zt*(4.*eos041))) )) )
503 drho_dt_dp = (d2r_p0 + zp*(d2r_p1 + zp*d2r_p2))
504
506
507!> Calculate the partial derivatives of specific volume with temperature and salinity
508!! using the density polynomial fit EOS from Roquet et al. (2015).
509elemental subroutine calculate_specvol_derivs_elem_roquet_rho(this, T, S, pressure, dSV_dT, dSV_dS)
510 class(roquet_rho_eos), intent(in) :: this !< This EOS
511 real, intent(in) :: t !< Conservative temperature [degC]
512 real, intent(in) :: s !< Absolute salinity [g kg-1]
513 real, intent(in) :: pressure !< Pressure [Pa]
514 real, intent(inout) :: dsv_dt !< The partial derivative of specific volume with
515 !! potential temperature [m3 kg-1 degC-1]
516 real, intent(inout) :: dsv_ds !< The partial derivative of specific volume with
517 !! salinity [m3 kg-1 ppt-1]
518 ! Local variables
519 real :: rho ! In situ density [kg m-3]
520 real :: drho_dt ! Derivative of density with temperature [kg m-3 degC-1]
521 real :: drho_ds ! Derivative of density with salinity [kg m-3 ppt-1]
522
523 call this%calculate_density_derivs_elem(t, s, pressure, drho_dt, drho_ds)
524 rho = this%density_elem(t, s, pressure)
525 dsv_dt = -drho_dt/(rho**2)
526 dsv_ds = -drho_ds/(rho**2)
527
529
530!> Compute the in situ density of sea water (rho in [kg m-3]) and the compressibility
531!! (drho/dp = C_sound^-2, stored as drho_dp [s2 m-2]) from absolute salinity (sal [g kg-1]),
532!! conservative temperature (T [degC]), and pressure [Pa], using the density polynomial
533!! fit EOS from Roquet et al. (2015).
534elemental subroutine calculate_compress_elem_roquet_rho(this, T, S, pressure, rho, drho_dp)
535 class(roquet_rho_eos), intent(in) :: this !< This EOS
536 real, intent(in) :: t !< Conservative temperature [degC]
537 real, intent(in) :: s !< Absolute salinity [g kg-1]
538 real, intent(in) :: pressure !< Pressure [Pa]
539 real, intent(out) :: rho !< In situ density [kg m-3]
540 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
541 !! (also the inverse of the square of sound speed)
542 !! [s2 m-2]
543 ! Local variables
544 real :: zp ! Pressure [Pa]
545 real :: zt ! Conservative temperature [degC]
546 real :: zs ! The square root of absolute salinity with an offset normalized
547 ! by an assumed salinity range [nondim]
548 real :: drho00p_dp ! Derivative of the pressure-dependent reference density profile with pressure [kg m-3 Pa-1]
549 real :: drhots_dp ! Derivative of the density anomaly from the reference profile with pressure [kg m-3 Pa-1]
550 real :: rho00p ! The pressure-dependent (but temperature and salinity independent) reference
551 ! density profile [kg m-3]
552 real :: rhots ! Density anomaly from the reference profile [kg m-3]
553 real :: rhots0 ! A contribution to density from temperature and salinity anomalies at the
554 ! surface pressure [kg m-3]
555 real :: rhots1 ! A density contribution proportional to pressure [kg m-3 Pa-1]
556 real :: rhots2 ! A density contribution proportional to pressure**2 [kg m-3 Pa-2]
557 real :: rhots3 ! A density contribution proportional to pressure**3 [kg m-3 Pa-3]
558 real :: rho0s0 ! Salinity dependent density at the surface pressure and zero temperature [kg m-3]
559
560 ! The following algorithm was published by Roquet et al. (2015), intended for use with NEMO.
561 ! Conversions to the units used here.
562 zt = t
563 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
564 zp = pressure
565
566 ! The next two lines should be used if it is necessary to convert potential temperature and
567 ! practical salinity to conservative temperature and absolute salinity.
568 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
569 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
570
571 rhots3 = eos003 + (zs*eos103 + zt*eos013)
572 rhots2 = eos002 + (zs*(eos102 + zs*eos202) &
573 + zt*(eos012 + (zs*eos112 + zt*eos022)) )
574 rhots1 = eos001 + (zs*(eos101 + zs*(eos201 + zs*(eos301 + zs*eos401))) &
575 + zt*(eos011 + (zs*(eos111 + zs*(eos211 + zs*eos311)) &
576 + zt*(eos021 + (zs*(eos121 + zs*eos221) &
577 + zt*(eos031 + (zs*eos131 + zt*eos041)) )) )) )
578
579 rhots0 = zt*(eos010 &
580 + (zs*(eos110 + zs*(eos210 + zs*(eos310 + zs*(eos410 + zs*eos510)))) &
581 + zt*(eos020 + (zs*(eos120 + zs*(eos220 + zs*(eos320 + zs*eos420))) &
582 + zt*(eos030 + (zs*(eos130 + zs*(eos230 + zs*eos330)) &
583 + zt*(eos040 + (zs*(eos140 + zs*eos240) &
584 + zt*(eos050 + (zs*eos150 + zt*eos060)) )) )) )) ) )
585
586 rho0s0 = eos000 + zs*(eos100 + zs*(eos200 + zs*(eos300 + zs*(eos400 + zs*(eos500 + zs*eos600)))))
587
588 rho00p = zp*(r00 + zp*(r01 + zp*(r02 + zp*(r03 + zp*(r04 + zp*r05)))))
589
590 rhots = (rhots0 + rho0s0) + zp*(rhots1 + zp*(rhots2 + zp*rhots3))
591 rho = rhots + rho00p ! In situ density [kg m-3]
592
593 drho00p_dp = r00 + zp*(2.*r01 + zp*(3.*r02 + zp*(4.*r03 + zp*(5.*r04 + zp*(6.*r05)))))
594 drhots_dp = rhots1 + zp*(2.*rhots2 + zp*(3.*rhots3))
595 drho_dp = drhots_dp + drho00p_dp ! Compressibility [s2 m-2]
596
598
599!> Return the range of temperatures, salinities and pressures for which the Roquet et al. (2015)
600!! expression for in situ density has been fitted to observations. Care should be taken when
601!! applying this equation of state outside of its fit range.
602subroutine eos_fit_range_roquet_rho(this, T_min, T_max, S_min, S_max, p_min, p_max)
603 class(roquet_rho_eos), intent(in) :: this !< This EOS
604 real, optional, intent(out) :: T_min !< The minimum conservative temperature over which this EoS is fitted [degC]
605 real, optional, intent(out) :: T_max !< The maximum conservative temperature over which this EoS is fitted [degC]
606 real, optional, intent(out) :: S_min !< The minimum absolute salinity over which this EoS is fitted [g kg-1]
607 real, optional, intent(out) :: S_max !< The maximum absolute salinity over which this EoS is fitted [g kg-1]
608 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
609 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
610
611 if (present(t_min)) t_min = -6.0
612 if (present(t_max)) t_max = 40.0
613 if (present(s_min)) s_min = 0.0
614 if (present(s_max)) s_max = 42.0
615 if (present(p_min)) p_min = 0.0
616 if (present(p_max)) p_max = 1.0e8
617
618end subroutine eos_fit_range_roquet_rho
619
620!> Calculate the in-situ density for 1D arraya inputs and outputs.
621subroutine calculate_density_array_roquet_rho(this, T, S, pressure, rho, start, npts, rho_ref)
622 class(roquet_rho_eos), intent(in) :: this !< This EOS
623 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
624 real, dimension(:), intent(in) :: S !< Salinity [PSU]
625 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
626 real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]
627 integer, intent(in) :: start !< The starting index for calculations
628 integer, intent(in) :: npts !< The number of values to calculate
629 real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]
630
631 ! Local variables
632 integer :: j
633
634 if (present(rho_ref)) then
635 do j = start, start+npts-1
636 rho(j) = density_anomaly_elem_roquet_rho(this, t(j), s(j), pressure(j), rho_ref)
637 enddo
638 else
639 do j = start, start+npts-1
640 rho(j) = density_elem_roquet_rho(this, t(j), s(j), pressure(j))
641 enddo
642 endif
643
645
646!> Calculate the in-situ specific volume for 1D array inputs and outputs.
647subroutine calculate_spec_vol_array_roquet_rho(this, T, S, pressure, specvol, start, npts, spv_ref)
648 class(roquet_rho_eos), intent(in) :: this !< This EOS
649 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
650 real, dimension(:), intent(in) :: S !< Salinity [PSU]
651 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
652 real, dimension(:), intent(out) :: specvol !< In situ specific volume [m3 kg-1]
653 integer, intent(in) :: start !< The starting index for calculations
654 integer, intent(in) :: npts !< The number of values to calculate
655 real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
656
657 ! Local variables
658 integer :: j
659
660 if (present(spv_ref)) then
661 do j = start, start+npts-1
662 specvol(j) = spec_vol_anomaly_elem_roquet_rho(this, t(j), s(j), pressure(j), spv_ref)
663 enddo
664 else
665 do j = start, start+npts-1
666 specvol(j) = spec_vol_elem_roquet_rho(this, t(j), s(j), pressure(j) )
667 enddo
668 endif
669
671
672!> \namespace mom_eos_Roquet_rho
673!!
674!! \section section_EOS_Roquet_rho Roquet_rho equation of state
675!!
676!! Fabien Roquet and colleagues developed this equation of state using a simple polynomial fit
677!! to the TEOS-10 equation of state, for efficiency when used in the NEMO ocean model. Fabien
678!! Roquet also graciously provided the MOM6 team with the original code implementing this
679!! equation of state, although it has since been modified and extended to have capabilities
680!! mirroring those available with other equations of state in MOM6. This particular equation
681!! of state is a balance between an accuracy that matches the TEOS-10 density to better than
682!! observational uncertainty with a polynomial form that can be evaluated quickly despite having
683!! 52 terms.
684!!
685!! \subsection section_EOS_Roquet_rho_references References
686!!
687!! Roquet, F., Madec, G., McDougall, T. J., and Barker, P. M., 2015:
688!! Accurate polynomial expressions for the density and specific volume
689!! of seawater using the TEOS-10 standard. Ocean Modelling, 90:29-43.
690
691end module mom_eos_roquet_rho