MOM_EOS_Roquet_SpV.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 for specific volume (SpV) using the expressions of Roquet et al. 2015
7
9
10implicit none ; private
11
12public roquet_spv_eos
13
14real, parameter :: pa2kb = 1.e-8 !< Conversion factor between Pa and kbar [kbar Pa-1]
15!>@{ Parameters in the Roquet specific volume polynomial equation of state
16real, parameter :: rdeltas = 24. ! 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! The following are the coefficients of the fit to the reference density profile (rho00p) as a function of
20! pressure (P), with a contribution R0c * P**(c+1). The nomenclature follows Roquet.
21real, parameter :: v00 = -4.4015007269e-05*pa2kb ! SpV00p P coef. [m3 kg-1 Pa-1]
22real, parameter :: v01 = 6.9232335784e-06*pa2kb**2 ! SpV00p P**2 coef. [m3 kg-1 Pa-2]
23real, parameter :: v02 = -7.5004675975e-07*pa2kb**3 ! SpV00p P**3 coef. [m3 kg-1 Pa-3]
24real, parameter :: v03 = 1.7009109288e-08*pa2kb**4 ! SpV00p P**4 coef. [m3 kg-1 Pa-4]
25real, parameter :: v04 = -1.6884162004e-08*pa2kb**5 ! SpV00p P**5 coef. [m3 kg-1 Pa-5]
26real, parameter :: v05 = 1.9613503930e-09*pa2kb**6 ! SpV00p P**6 coef. [m3 kg-1 Pa-6]
27
28! The following terms are contributions to specific volume (SpV) as a function of the square root of
29! normalized absolute salinity with an offset (zs), temperature (T) and pressure (P), with a contribution
30! SPVabc * zs**a * T**b * P**c. The numbers here are copied directly from Roquet et al. (2015), but
31! the expressions here do not use the same nondimensionalization for pressure or temperature as they do.
32real, parameter :: spv000 = 1.0772899069e-03 ! Constant SpV contribution [m3 kg-1]
33real, parameter :: spv100 = -3.1263658781e-04 ! SpV zs coef. [m3 kg-1]
34real, parameter :: spv200 = 6.7615860683e-04 ! SpV zs**2 coef. [m3 kg-1]
35real, parameter :: spv300 = -8.6127884515e-04 ! SpV zs**3 coef. [m3 kg-1]
36real, parameter :: spv400 = 5.9010812596e-04 ! SpV zs**4 coef. [m3 kg-1]
37real, parameter :: spv500 = -2.1503943538e-04 ! SpV zs**5 coef. [m3 kg-1]
38real, parameter :: spv600 = 3.2678954455e-05 ! SpV zs**6 coef. [m3 kg-1]
39real, parameter :: spv010 = -1.4949652640e-05*i_ts ! SpV T coef. [m3 kg-1 degC-1]
40real, parameter :: spv110 = 3.1866349188e-05*i_ts ! SpV zs * T coef. [m3 kg-1 degC-1]
41real, parameter :: spv210 = -3.8070687610e-05*i_ts ! SpV zs**2 * T coef. [m3 kg-1 degC-1]
42real, parameter :: spv310 = 2.9818473563e-05*i_ts ! SpV zs**3 * T coef. [m3 kg-1 degC-1]
43real, parameter :: spv410 = -1.0011321965e-05*i_ts ! SpV zs**4 * T coef. [m3 kg-1 degC-1]
44real, parameter :: spv510 = 1.0751931163e-06*i_ts ! SpV zs**5 * T coef. [m3 kg-1 degC-1]
45real, parameter :: spv020 = 2.7546851539e-05*i_ts**2 ! SpV T**2 coef. [m3 kg-1 degC-2]
46real, parameter :: spv120 = -3.6597334199e-05*i_ts**2 ! SpV zs * T**2 coef. [m3 kg-1 degC-2]
47real, parameter :: spv220 = 3.4489154625e-05*i_ts**2 ! SpV zs**2 * T**2 coef. [m3 kg-1 degC-2]
48real, parameter :: spv320 = -1.7663254122e-05*i_ts**2 ! SpV zs**3 * T**2 coef. [m3 kg-1 degC-2]
49real, parameter :: spv420 = 3.5965131935e-06*i_ts**2 ! SpV zs**4 * T**2 coef. [m3 kg-1 degC-2]
50real, parameter :: spv030 = -1.6506828994e-05*i_ts**3 ! SpV T**3 coef. [m3 kg-1 degC-3]
51real, parameter :: spv130 = 2.4412359055e-05*i_ts**3 ! SpV zs * T**3 coef. [m3 kg-1 degC-3]
52real, parameter :: spv230 = -1.4606740723e-05*i_ts**3 ! SpV zs**2 * T**3 coef. [m3 kg-1 degC-3]
53real, parameter :: spv330 = 2.3293406656e-06*i_ts**3 ! SpV zs**3 * T**3 coef. [m3 kg-1 degC-3]
54real, parameter :: spv040 = 6.7896174634e-06*i_ts**4 ! SpV T**4 coef. [m3 kg-1 degC-4]
55real, parameter :: spv140 = -8.7951832993e-06*i_ts**4 ! SpV zs * T**4 coef. [m3 kg-1 degC-4]
56real, parameter :: spv240 = 4.4249040774e-06*i_ts**4 ! SpV zs**2 * T**4 coef. [m3 kg-1 degC-4]
57real, parameter :: spv050 = -7.2535743349e-07*i_ts**5 ! SpV T**5 coef. [m3 kg-1 degC-5]
58real, parameter :: spv150 = -3.4680559205e-07*i_ts**5 ! SpV zs * T**5 coef. [m3 kg-1 degC-5]
59real, parameter :: spv060 = 1.9041365570e-07*i_ts**6 ! SpV T**6 coef. [m3 kg-1 degC-6]
60real, parameter :: spv001 = -1.6889436589e-05*pa2kb ! SpV P coef. [m3 kg-1 Pa-1]
61real, parameter :: spv101 = 2.1106556158e-05*pa2kb ! SpV zs * P coef. [m3 kg-1 Pa-1]
62real, parameter :: spv201 = -2.1322804368e-05*pa2kb ! SpV zs**2 * P coef. [m3 kg-1 Pa-1]
63real, parameter :: spv301 = 1.7347655458e-05*pa2kb ! SpV zs**3 * P coef. [m3 kg-1 Pa-1]
64real, parameter :: spv401 = -4.3209400767e-06*pa2kb ! SpV zs**4 * P coef. [m3 kg-1 Pa-1]
65real, parameter :: spv011 = 1.5355844621e-05*(i_ts*pa2kb) ! SpV T * P coef. [m3 kg-1 degC-1 Pa-1]
66real, parameter :: spv111 = 2.0914122241e-06*(i_ts*pa2kb) ! SpV zs * T * P coef. [m3 kg-1 degC-1 Pa-1]
67real, parameter :: spv211 = -5.7751479725e-06*(i_ts*pa2kb) ! SpV zs**2 * T * P coef. [m3 kg-1 degC-1 Pa-1]
68real, parameter :: spv311 = 1.0767234341e-06*(i_ts*pa2kb) ! SpV zs**3 * T * P coef. [m3 kg-1 degC-1 Pa-1]
69real, parameter :: spv021 = -9.6659393016e-06*(i_ts**2*pa2kb) ! SpV T**2 * P coef. [m3 kg-1 degC-2 Pa-1]
70real, parameter :: spv121 = -7.0686982208e-07*(i_ts**2*pa2kb) ! SpV zs * T**2 * P coef. [m3 kg-1 degC-2 Pa-1]
71real, parameter :: spv221 = 1.4488066593e-06*(i_ts**2*pa2kb) ! SpV zs**2 * T**2 * P coef. [m3 kg-1 degC-2 Pa-1]
72real, parameter :: spv031 = 3.1134283336e-06*(i_ts**3*pa2kb) ! SpV T**3 * P coef. [m3 kg-1 degC-3 Pa-1]
73real, parameter :: spv131 = 7.9562529879e-08*(i_ts**3*pa2kb) ! SpV zs * T**3 * P coef. [m3 kg-1 degC-3 Pa-1]
74real, parameter :: spv041 = -5.6590253863e-07*(i_ts**4*pa2kb) ! SpV T**4 * P coef. [m3 kg-1 degC-4 Pa-1]
75real, parameter :: spv002 = 1.0500241168e-06*pa2kb**2 ! SpV P**2 coef. [m3 kg-1 Pa-2]
76real, parameter :: spv102 = 1.9600661704e-06*pa2kb**2 ! SpV zs * P**2 coef. [m3 kg-1 Pa-2]
77real, parameter :: spv202 = -2.1666693382e-06*pa2kb**2 ! SpV zs**2 * P**2 coef. [m3 kg-1 Pa-2]
78real, parameter :: spv012 = -3.8541359685e-06*(i_ts*pa2kb**2) ! SpV T * P**2 coef. [m3 kg-1 degC-1 Pa-2]
79real, parameter :: spv112 = 1.0157632247e-06*(i_ts*pa2kb**2) ! SpV zs * T * P**2 coef. [m3 kg-1 degC-1 Pa-2]
80real, parameter :: spv022 = 1.7178343158e-06*(i_ts**2*pa2kb**2) ! SpV T**2 * P**2 coef. [m3 kg-1 degC-2 Pa-2]
81real, parameter :: spv003 = -4.1503454190e-07*pa2kb**3 ! SpV P**3 coef. [m3 kg-1 Pa-3]
82real, parameter :: spv103 = 3.5627020989e-07*pa2kb**3 ! SpV zs * P**3 coef. [m3 kg-1 Pa-3]
83real, parameter :: spv013 = -1.1293871415e-07*(i_ts*pa2kb**3) ! SpV T * P**3 coef. [m3 kg-1 degC-1 Pa-3]
84
85real, parameter :: alp000 = spv010 ! Constant in the dSpV_dT fit [m3 kg-1 degC-1]
86real, parameter :: alp100 = spv110 ! dSpV_dT fit zs coef. [m3 kg-1 degC-1]
87real, parameter :: alp200 = spv210 ! dSpV_dT fit zs**2 coef. [m3 kg-1 degC-1]
88real, parameter :: alp300 = spv310 ! dSpV_dT fit zs**3 coef. [m3 kg-1 degC-1]
89real, parameter :: alp400 = spv410 ! dSpV_dT fit zs**4 coef. [m3 kg-1 degC-1]
90real, parameter :: alp500 = spv510 ! dSpV_dT fit zs**5 coef. [m3 kg-1 degC-1]
91real, parameter :: alp010 = 2.*spv020 ! dSpV_dT fit T coef. [m3 kg-1 degC-2]
92real, parameter :: alp110 = 2.*spv120 ! dSpV_dT fit zs * T coef. [m3 kg-1 degC-2]
93real, parameter :: alp210 = 2.*spv220 ! dSpV_dT fit zs**2 * T coef. [m3 kg-1 degC-2]
94real, parameter :: alp310 = 2.*spv320 ! dSpV_dT fit zs**3 * T coef. [m3 kg-1 degC-2]
95real, parameter :: alp410 = 2.*spv420 ! dSpV_dT fit zs**4 * T coef. [m3 kg-1 degC-2]
96real, parameter :: alp020 = 3.*spv030 ! dSpV_dT fit T**2 coef. [m3 kg-1 degC-3]
97real, parameter :: alp120 = 3.*spv130 ! dSpV_dT fit zs * T**2 coef. [m3 kg-1 degC-3]
98real, parameter :: alp220 = 3.*spv230 ! dSpV_dT fit zs**2 * T**2 coef. [m3 kg-1 degC-3]
99real, parameter :: alp320 = 3.*spv330 ! dSpV_dT fit zs**3 * T**2 coef. [m3 kg-1 degC-3]
100real, parameter :: alp030 = 4.*spv040 ! dSpV_dT fit T**3 coef. [m3 kg-1 degC-4]
101real, parameter :: alp130 = 4.*spv140 ! dSpV_dT fit zs * T**3 coef. [m3 kg-1 degC-4]
102real, parameter :: alp230 = 4.*spv240 ! dSpV_dT fit zs**2 * T**3 coef. [m3 kg-1 degC-4]
103real, parameter :: alp040 = 5.*spv050 ! dSpV_dT fit T**4 coef. [m3 kg-1 degC-5]
104real, parameter :: alp140 = 5.*spv150 ! dSpV_dT fit zs* * T**4 coef. [m3 kg-1 degC-5]
105real, parameter :: alp050 = 6.*spv060 ! dSpV_dT fit T**5 coef. [m3 kg-1 degC-6]
106real, parameter :: alp001 = spv011 ! dSpV_dT fit P coef. [m3 kg-1 degC-1 Pa-1]
107real, parameter :: alp101 = spv111 ! dSpV_dT fit zs * P coef. [m3 kg-1 degC-1 Pa-1]
108real, parameter :: alp201 = spv211 ! dSpV_dT fit zs**2 * P coef. [m3 kg-1 degC-1 Pa-1]
109real, parameter :: alp301 = spv311 ! dSpV_dT fit zs**3 * P coef. [m3 kg-1 degC-1 Pa-1]
110real, parameter :: alp011 = 2.*spv021 ! dSpV_dT fit T * P coef. [m3 kg-1 degC-2 Pa-1]
111real, parameter :: alp111 = 2.*spv121 ! dSpV_dT fit zs * T * P coef. [m3 kg-1 degC-2 Pa-1]
112real, parameter :: alp211 = 2.*spv221 ! dSpV_dT fit zs**2 * T * P coef. [m3 kg-1 degC-2 Pa-1]
113real, parameter :: alp021 = 3.*spv031 ! dSpV_dT fit T**2 * P coef. [m3 kg-1 degC-3 Pa-1]
114real, parameter :: alp121 = 3.*spv131 ! dSpV_dT fit zs * T**2 * P coef. [m3 kg-1 degC-3 Pa-1]
115real, parameter :: alp031 = 4.*spv041 ! dSpV_dT fit T**3 * P coef. [m3 kg-1 degC-4 Pa-1]
116real, parameter :: alp002 = spv012 ! dSpV_dT fit P**2 coef. [m3 kg-1 degC-1 Pa-2]
117real, parameter :: alp102 = spv112 ! dSpV_dT fit zs * P**2 coef. [m3 kg-1 degC-1 Pa-2]
118real, parameter :: alp012 = 2.*spv022 ! dSpV_dT fit T * P**2 coef. [m3 kg-1 degC-2 Pa-2]
119real, parameter :: alp003 = spv013 ! dSpV_dT fit P**3 coef. [m3 kg-1 degC-1 Pa-3]
120
121real, parameter :: bet000 = 0.5*spv100*r1_s0 ! Constant in the dSpV_dS fit [m3 kg-1 ppt-1]
122real, parameter :: bet100 = spv200*r1_s0 ! dSpV_dS fit zs coef. [m3 kg-1 ppt-1]
123real, parameter :: bet200 = 1.5*spv300*r1_s0 ! dSpV_dS fit zs**2 coef. [m3 kg-1 ppt-1]
124real, parameter :: bet300 = 2.0*spv400*r1_s0 ! dSpV_dS fit zs**3 coef. [m3 kg-1 ppt-1]
125real, parameter :: bet400 = 2.5*spv500*r1_s0 ! dSpV_dS fit zs**4 coef. [m3 kg-1 ppt-1]
126real, parameter :: bet500 = 3.0*spv600*r1_s0 ! dSpV_dS fit zs**5 coef. [m3 kg-1 ppt-1]
127real, parameter :: bet010 = 0.5*spv110*r1_s0 ! dSpV_dS fit T coef. [m3 kg-1 ppt-1 degC-1]
128real, parameter :: bet110 = spv210*r1_s0 ! dSpV_dS fit zs * T coef. [m3 kg-1 ppt-1 degC-1]
129real, parameter :: bet210 = 1.5*spv310*r1_s0 ! dSpV_dS fit zs**2 * T coef. [m3 kg-1 ppt-1 degC-1]
130real, parameter :: bet310 = 2.0*spv410*r1_s0 ! dSpV_dS fit zs**3 * T coef. [m3 kg-1 ppt-1 degC-1]
131real, parameter :: bet410 = 2.5*spv510*r1_s0 ! dSpV_dS fit zs**4 * T coef. [m3 kg-1 ppt-1 degC-1]
132real, parameter :: bet020 = 0.5*spv120*r1_s0 ! dSpV_dS fit T**2 coef. [m3 kg-1 ppt-1 degC-2]
133real, parameter :: bet120 = spv220*r1_s0 ! dSpV_dS fit zs * T**2 coef. [m3 kg-1 ppt-1 degC-2]
134real, parameter :: bet220 = 1.5*spv320*r1_s0 ! dSpV_dS fit zs**2 * T**2 coef. [m3 kg-1 ppt-1 degC-2]
135real, parameter :: bet320 = 2.0*spv420*r1_s0 ! dSpV_dS fit zs**3 * T**2 coef. [m3 kg-1 ppt-1 degC-2]
136real, parameter :: bet030 = 0.5*spv130*r1_s0 ! dSpV_dS fit T**3 coef. [m3 kg-1 ppt-1 degC-3]
137real, parameter :: bet130 = spv230*r1_s0 ! dSpV_dS fit zs * T**3 coef. [m3 kg-1 ppt-1 degC-3]
138real, parameter :: bet230 = 1.5*spv330*r1_s0 ! dSpV_dS fit zs**2 * T**3 coef. [m3 kg-1 ppt-1 degC-3]
139real, parameter :: bet040 = 0.5*spv140*r1_s0 ! dSpV_dS fit T**4 coef. [m3 kg-1 ppt-1 degC-4]
140real, parameter :: bet140 = spv240*r1_s0 ! dSpV_dS fit zs * T**4 coef. [m3 kg-1 ppt-1 degC-4]
141real, parameter :: bet050 = 0.5*spv150*r1_s0 ! dSpV_dS fit T**5 coef. [m3 kg-1 ppt-1 degC-5]
142real, parameter :: bet001 = 0.5*spv101*r1_s0 ! dSpV_dS fit P coef. [m3 kg-1 ppt-1 Pa-1]
143real, parameter :: bet101 = spv201*r1_s0 ! dSpV_dS fit zs * P coef. [m3 kg-1 ppt-1 Pa-1]
144real, parameter :: bet201 = 1.5*spv301*r1_s0 ! dSpV_dS fit zs**2 * P coef. [m3 kg-1 ppt-1 Pa-1]
145real, parameter :: bet301 = 2.0*spv401*r1_s0 ! dSpV_dS fit zs**3 * P coef. [m3 kg-1 ppt-1 Pa-1]
146real, parameter :: bet011 = 0.5*spv111*r1_s0 ! dSpV_dS fit T * P coef. [m3 kg-1 ppt-1 degC-1 Pa-1]
147real, parameter :: bet111 = spv211*r1_s0 ! dSpV_dS fit zs * T * P coef. [m3 kg-1 ppt-1 degC-1 Pa-1]
148real, parameter :: bet211 = 1.5*spv311*r1_s0 ! dSpV_dS fit zs**2 * T * P coef. [m3 kg-1 ppt-1 degC-1 Pa-1]
149real, parameter :: bet021 = 0.5*spv121*r1_s0 ! dSpV_dS fit T**2 * P coef. [m3 kg-1 ppt-1 degC-2 Pa-1]
150real, parameter :: bet121 = spv221*r1_s0 ! dSpV_dS fit zs * T**2 * P coef. [m3 kg-1 ppt-1 degC-2 Pa-1]
151real, parameter :: bet031 = 0.5*spv131*r1_s0 ! dSpV_dS fit T**3 * P coef. [m3 kg-1 ppt-1 degC-3 Pa-1]
152real, parameter :: bet002 = 0.5*spv102*r1_s0 ! dSpV_dS fit P**2 coef. [m3 kg-1 ppt-1 Pa-2]
153real, parameter :: bet102 = spv202*r1_s0 ! dSpV_dS fit zs * P**2 coef. [m3 kg-1 ppt-1 Pa-2]
154real, parameter :: bet012 = 0.5*spv112*r1_s0 ! dSpV_dS fit T * P**2 coef. [m3 kg-1 ppt-1 degC-1 Pa-2]
155real, parameter :: bet003 = 0.5*spv103*r1_s0 ! dSpV_dS fit P**3 coef. [m3 kg-1 ppt-1 Pa-3]
156!>@}
157
158!> The EOS_base implementation of the Roquet et al., 2015, equation of state
159type, extends (eos_base) :: roquet_spv_eos
160
161contains
162 !> Implementation of the in-situ density as an elemental function [kg m-3]
163 procedure :: density_elem => density_elem_roquet_spv
164 !> Implementation of the in-situ density anomaly as an elemental function [kg m-3]
165 procedure :: density_anomaly_elem => density_anomaly_elem_roquet_spv
166 !> Implementation of the in-situ specific volume as an elemental function [m3 kg-1]
167 procedure :: spec_vol_elem => spec_vol_elem_roquet_spv
168 !> Implementation of the in-situ specific volume anomaly as an elemental function [m3 kg-1]
169 procedure :: spec_vol_anomaly_elem => spec_vol_anomaly_elem_roquet_spv
170 !> Implementation of the calculation of derivatives of density
171 procedure :: calculate_density_derivs_elem => calculate_density_derivs_elem_roquet_spv
172 !> Implementation of the calculation of second derivatives of density
173 procedure :: calculate_density_second_derivs_elem => calculate_density_second_derivs_elem_roquet_spv
174 !> Implementation of the calculation of derivatives of specific volume
175 procedure :: calculate_specvol_derivs_elem => calculate_specvol_derivs_elem_roquet_spv
176 !> Implementation of the calculation of compressibility
177 procedure :: calculate_compress_elem => calculate_compress_elem_roquet_spv
178 !> Implementation of the range query function
179 procedure :: eos_fit_range => eos_fit_range_roquet_spv
180
181 !> Local implementation of generic calculate_density_array for efficiency
182 procedure :: calculate_density_array => calculate_density_array_roquet_spv
183 !> Local implementation of generic calculate_spec_vol_array for efficiency
184 procedure :: calculate_spec_vol_array => calculate_spec_vol_array_roquet_spv
185
186end type roquet_spv_eos
187
188contains
189
190!> Roquet et al. in situ specific volume of sea water [m3 kg-1]
191!!
192!! This is an elemental function that can be applied to any combination of scalar and array inputs.
193real elemental function spec_vol_elem_roquet_spv(this, t, s, pressure)
194 class(roquet_spv_eos), intent(in) :: this !< This EOS
195 real, intent(in) :: t !< Conservative temperature [degC]
196 real, intent(in) :: s !< Absolute salinity [g kg-1]
197 real, intent(in) :: pressure !< pressure [Pa]
198
199 ! Local variables
200 real :: zp ! Pressure [Pa]
201 real :: zt ! Conservative temperature [degC]
202 real :: zs ! The square root of absolute salinity with an offset normalized
203 ! by an assumed salinity range [nondim]
204 real :: sv_00p ! A pressure-dependent but temperature and salinity independent contribution to
205 ! specific volume at the reference temperature and salinity [m3 kg-1]
206 real :: sv_ts ! Specific volume without a pressure-dependent contribution [m3 kg-1]
207 real :: sv_ts0 ! A contribution to specific volume from temperature and salinity anomalies at
208 ! the surface pressure [m3 kg-1]
209 real :: sv_ts1 ! A temperature and salinity dependent specific volume contribution that is
210 ! proportional to pressure [m3 kg-1 Pa-1]
211 real :: sv_ts2 ! A temperature and salinity dependent specific volume contribution that is
212 ! proportional to pressure**2 [m3 kg-1 Pa-2]
213 real :: sv_ts3 ! A temperature and salinity dependent specific volume contribution that is
214 ! proportional to pressure**3 [m3 kg-1 Pa-3]
215 real :: sv_0s0 ! Salinity dependent specific volume at the surface pressure and zero temperature [m3 kg-1]
216
217 ! The following algorithm was published by Roquet et al. (2015), intended for use in non-Boussinesq ocean models.
218
219 ! Conversions to the units used here.
220 zt = t
221 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
222 zp = pressure
223
224 ! The next two lines should be used if it is necessary to convert potential temperature and
225 ! practical salinity to conservative temperature and absolute salinity.
226 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
227 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
228
229 sv_ts3 = spv003 + (zs*spv103 + zt*spv013)
230 sv_ts2 = spv002 + (zs*(spv102 + zs*spv202) &
231 + zt*(spv012 + (zs*spv112 + zt*spv022)) )
232 sv_ts1 = spv001 + (zs*(spv101 + zs*(spv201 + zs*(spv301 + zs*spv401))) &
233 + zt*(spv011 + (zs*(spv111 + zs*(spv211 + zs*spv311)) &
234 + zt*(spv021 + (zs*(spv121 + zs*spv221) &
235 + zt*(spv031 + (zs*spv131 + zt*spv041)) )) )) )
236 sv_ts0 = zt*(spv010 &
237 + (zs*(spv110 + zs*(spv210 + zs*(spv310 + zs*(spv410 + zs*spv510)))) &
238 + zt*(spv020 + (zs*(spv120 + zs*(spv220 + zs*(spv320 + zs*spv420))) &
239 + zt*(spv030 + (zs*(spv130 + zs*(spv230 + zs*spv330)) &
240 + zt*(spv040 + (zs*(spv140 + zs*spv240) &
241 + zt*(spv050 + (zs*spv150 + zt*spv060)) )) )) )) ) )
242
243 sv_0s0 = spv000 + zs*(spv100 + zs*(spv200 + zs*(spv300 + zs*(spv400 + zs*(spv500 + zs*spv600)))))
244
245 sv_00p = zp*(v00 + zp*(v01 + zp*(v02 + zp*(v03 + zp*(v04 + zp*v05)))))
246
247 sv_ts = (sv_ts0 + sv_0s0) + zp*(sv_ts1 + zp*(sv_ts2 + zp*sv_ts3))
248 spec_vol_elem_roquet_spv = sv_ts + sv_00p ! In situ specific volume [m3 kg-1]
249
250end function spec_vol_elem_roquet_spv
251
252!> Roquet et al. in situ specific volume anomaly of sea water [m3 kg-1]
253!!
254!! This is an elemental function that can be applied to any combination of scalar and array inputs.
255real elemental function spec_vol_anomaly_elem_roquet_spv(this, t, s, pressure, spv_ref)
256 class(roquet_spv_eos), intent(in) :: this !< This EOS
257 real, intent(in) :: t !< Conservative temperature [degC]
258 real, intent(in) :: s !< Absolute salinity [g kg-1]
259 real, intent(in) :: pressure !< pressure [Pa]
260 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
261
262 ! Local variables
263 real :: zp ! Pressure [Pa]
264 real :: zt ! Conservative temperature [degC]
265 real :: zs ! The square root of absolute salinity with an offset normalized
266 ! by an assumed salinity range [nondim]
267 real :: sv_00p ! A pressure-dependent but temperature and salinity independent contribution to
268 ! specific volume at the reference temperature and salinity [m3 kg-1]
269 real :: sv_ts ! Specific volume without a pressure-dependent contribution [m3 kg-1]
270 real :: sv_ts0 ! A contribution to specific volume from temperature and salinity anomalies at
271 ! the surface pressure [m3 kg-1]
272 real :: sv_ts1 ! A temperature and salinity dependent specific volume contribution that is
273 ! proportional to pressure [m3 kg-1 Pa-1]
274 real :: sv_ts2 ! A temperature and salinity dependent specific volume contribution that is
275 ! proportional to pressure**2 [m3 kg-1 Pa-2]
276 real :: sv_ts3 ! A temperature and salinity dependent specific volume contribution that is
277 ! proportional to pressure**3 [m3 kg-1 Pa-3]
278 real :: sv_0s0 ! Salinity dependent specific volume at the surface pressure and zero temperature [m3 kg-1]
279
280 ! The following algorithm was published by Roquet et al. (2015), intended for use in non-Boussinesq ocean models.
281
282 ! Conversions to the units used here.
283 zt = t
284 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
285 zp = pressure
286
287 ! The next two lines should be used if it is necessary to convert potential temperature and
288 ! practical salinity to conservative temperature and absolute salinity.
289 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
290 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
291
292 sv_ts3 = spv003 + (zs*spv103 + zt*spv013)
293 sv_ts2 = spv002 + (zs*(spv102 + zs*spv202) &
294 + zt*(spv012 + (zs*spv112 + zt*spv022)) )
295 sv_ts1 = spv001 + (zs*(spv101 + zs*(spv201 + zs*(spv301 + zs*spv401))) &
296 + zt*(spv011 + (zs*(spv111 + zs*(spv211 + zs*spv311)) &
297 + zt*(spv021 + (zs*(spv121 + zs*spv221) &
298 + zt*(spv031 + (zs*spv131 + zt*spv041)) )) )) )
299 sv_ts0 = zt*(spv010 &
300 + (zs*(spv110 + zs*(spv210 + zs*(spv310 + zs*(spv410 + zs*spv510)))) &
301 + zt*(spv020 + (zs*(spv120 + zs*(spv220 + zs*(spv320 + zs*spv420))) &
302 + zt*(spv030 + (zs*(spv130 + zs*(spv230 + zs*spv330)) &
303 + zt*(spv040 + (zs*(spv140 + zs*spv240) &
304 + zt*(spv050 + (zs*spv150 + zt*spv060)) )) )) )) ) )
305
306 sv_0s0 = spv000 + zs*(spv100 + zs*(spv200 + zs*(spv300 + zs*(spv400 + zs*(spv500 + zs*spv600)))))
307
308 sv_00p = zp*(v00 + zp*(v01 + zp*(v02 + zp*(v03 + zp*(v04 + zp*v05)))))
309
310 sv_0s0 = sv_0s0 - spv_ref
311
312 sv_ts = (sv_ts0 + sv_0s0) + zp*(sv_ts1 + zp*(sv_ts2 + zp*sv_ts3))
313 spec_vol_anomaly_elem_roquet_spv = sv_ts + sv_00p ! In situ specific volume [m3 kg-1]
314
316
317!> Roquet in situ density [kg m-3]
318!!
319!! This is an elemental function that can be applied to any combination of scalar and array inputs.
320real elemental function density_elem_roquet_spv(this, t, s, pressure)
321 class(roquet_spv_eos), intent(in) :: this !< This EOS
322 real, intent(in) :: t !< Conservative temperature [degC]
323 real, intent(in) :: s !< Absolute salinity [g kg-1]
324 real, intent(in) :: pressure !< Pressure [Pa]
325
326 ! Local variables
327 real :: spv ! The specific volume [m3 kg-1]
328
329 spv = spec_vol_elem_roquet_spv(this, t, s, pressure)
330 density_elem_roquet_spv = 1.0 / spv ! In situ density [kg m-3]
331
332end function density_elem_roquet_spv
333
334!> Roquet in situ density anomaly [kg m-3]
335!!
336!! This is an elemental function that can be applied to any combination of scalar and array inputs.
337real elemental function density_anomaly_elem_roquet_spv(this, t, s, pressure, rho_ref)
338 class(roquet_spv_eos), intent(in) :: this !< This EOS
339 real, intent(in) :: t !< Conservative temperature [degC]
340 real, intent(in) :: s !< Absolute salinity [g kg-1]
341 real, intent(in) :: pressure !< Pressure [Pa]
342 real, intent(in) :: rho_ref !< A reference density [kg m-3]
343
344 ! Local variables
345 real :: spv ! The specific volume [m3 kg-1]
346
347 spv = spec_vol_anomaly_elem_roquet_spv(this, t, s, pressure, spv_ref=1.0/rho_ref)
348 density_anomaly_elem_roquet_spv = -rho_ref**2*spv / (rho_ref*spv + 1.0) ! In situ density [kg m-3]
349
351
352!> Return the partial derivatives of specific volume with temperature and salinity for 1-d array
353!! inputs and outputs, using the specific volume polynomial fit from Roquet et al. (2015).
354elemental subroutine calculate_specvol_derivs_elem_roquet_spv(this, T, S, pressure, dSV_dT, dSV_dS)
355 class(roquet_spv_eos), intent(in) :: this !< This EOS
356 real, intent(in) :: t !< Conservative temperature [degC]
357 real, intent(in) :: s !< Absolute salinity [g kg-1]
358 real, intent(in) :: pressure !< Pressure [Pa]
359 real, intent(inout) :: dsv_dt !< The partial derivative of specific volume with
360 !! conservative temperature [m3 kg-1 degC-1]
361 real, intent(inout) :: dsv_ds !< The partial derivative of specific volume with
362 !! absolute salinity [m3 kg-1 ppt-1]
363
364 real :: zp ! Pressure [Pa]
365 real :: zt ! Conservative temperature [degC]
366 real :: zs ! The square root of absolute salinity with an offset normalized
367 ! by an assumed salinity range [nondim]
368 real :: dsvdzt0 ! A contribution to the partial derivative of specific volume with temperature
369 ! from temperature anomalies at the surface pressure [m3 kg-1 degC-1]
370 real :: dsvdzt1 ! A contribution to the partial derivative of specific volume with temperature
371 ! that is proportional to pressure [m3 kg-1 degC-1 Pa-1]
372 real :: dsvdzt2 ! A contribution to the partial derivative of specific volume with temperature
373 ! that is proportional to pressure**2 [m3 kg-1 degC-1 Pa-2]
374 real :: dsvdzt3 ! A contribution to the partial derivative of specific volume with temperature
375 ! that is proportional to pressure**3 [m3 kg-1 degC-1 Pa-3]
376 real :: dsvdzs0 ! A contribution to the partial derivative of specific volume with
377 ! salinity [m3 kg-1 ppt-1] from temperature anomalies at the surface pressure
378 real :: dsvdzs1 ! A contribution to the partial derivative of specific volume with
379 ! salinity [m3 kg-1 ppt-1 Pa-1] proportional to pressure
380 real :: dsvdzs2 ! A contribution to the partial derivative of specific volume with
381 ! salinity [m3 kg-1 ppt-1 Pa-2] proportional to pressure**2
382 real :: dsvdzs3 ! A contribution to the partial derivative of specific volume with
383 ! salinity [m3 kg-1 ppt-1 Pa-3] proportional to pressure**3
384
385 ! Conversions to the units used here.
386 zt = t
387 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
388 zp = pressure
389
390 ! The next two lines should be used if it is necessary to convert potential temperature and
391 ! practical salinity to conservative temperature and absolute salinity.
392 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
393 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
394
395 ! Find the partial derivative of specific volume with temperature
396 dsvdzt3 = alp003
397 dsvdzt2 = alp002 + (zs*alp102 + zt*alp012)
398 dsvdzt1 = alp001 + (zs*(alp101 + zs*(alp201 + zs*alp301)) &
399 + zt*(alp011 + (zs*(alp111 + zs*alp211) &
400 + zt*(alp021 + (zs*alp121 + zt*alp031)) )) )
401 dsvdzt0 = alp000 + (zs*(alp100 + zs*(alp200 + zs*(alp300 + zs*(alp400 + zs*alp500)))) &
402 + zt*(alp010 + (zs*(alp110 + zs*(alp210 + zs*(alp310 + zs*alp410))) &
403 + zt*(alp020 + (zs*(alp120 + zs*(alp220 + zs*alp320)) &
404 + zt*(alp030 + (zt*(alp040 + (zs*alp140 + zt*alp050)) &
405 + zs*(alp130 + zs*alp230) )) )) )) )
406
407 dsv_dt = dsvdzt0 + zp*(dsvdzt1 + zp*(dsvdzt2 + zp*dsvdzt3))
408
409 ! Find the partial derivative of specific volume with salinity
410 dsvdzs3 = bet003
411 dsvdzs2 = bet002 + (zs*bet102 + zt*bet012)
412 dsvdzs1 = bet001 + (zs*(bet101 + zs*(bet201 + zs*bet301)) &
413 + zt*(bet011 + (zs*(bet111 + zs*bet211) &
414 + zt*(bet021 + (zs*bet121 + zt*bet031)) )) )
415 dsvdzs0 = bet000 + (zs*(bet100 + zs*(bet200 + zs*(bet300 + zs*(bet400 + zs*bet500)))) &
416 + zt*(bet010 + (zs*(bet110 + zs*(bet210 + zs*(bet310 + zs*bet410))) &
417 + zt*(bet020 + (zs*(bet120 + zs*(bet220 + zs*bet320)) &
418 + zt*(bet030 + (zt*(bet040 + (zs*bet140 + zt*bet050)) &
419 + zs*(bet130 + zs*bet230) )) )) )) )
420
421 ! The division by zs here is because zs = sqrt(S + S0), so dSV_dS = dzs_dS * dSV_dzs = (0.5 / zs) * dSV_dzs
422 dsv_ds = (dsvdzs0 + zp*(dsvdzs1 + zp*(dsvdzs2 + zp * dsvdzs3))) / zs
423
425
426!> Compute an array of derivatives of densities of sea water with temperature (drho_dT in [kg m-3 degC-1])
427!! and salinity (drho_dS in [kg m-3 ppt-1]) from absolute salinity (S [g kg-1]), conservative temperature
428!! (T [degC]) and pressure [Pa], using the specific volume polynomial fit from Roquet et al. (2015).
429elemental subroutine calculate_density_derivs_elem_roquet_spv(this, T, S, pressure, drho_dT, drho_dS)
430 class(roquet_spv_eos), intent(in) :: this !< This EOS
431 real, intent(in) :: t !< Conservative temperature [degC]
432 real, intent(in) :: s !< Absolute salinity [g kg-1]
433 real, intent(in) :: pressure !< pressure [Pa]
434 real, intent(out) :: drho_dt !< The partial derivative of density with
435 !! conservative temperature [kg m-3 degC-1]
436 real, intent(out) :: drho_ds !< The partial derivative of density with
437 !! absolute salinity [kg m-3 ppt-1]
438
439 ! Local variables
440 real :: dsv_dt ! The partial derivative of specific volume with
441 ! conservative temperature [m3 kg-1 degC-1]
442 real :: dsv_ds ! The partial derivative of specific volume with
443 ! absolute salinity [m3 kg-1 ppt-1]
444 real :: specvol ! The specific volume [m3 kg-1]
445 real :: rho ! The in situ density [kg m-3]
446
447 call this%calculate_specvol_derivs_elem(t, s, pressure, dsv_dt, dsv_ds)
448
449 specvol = this%spec_vol_elem(t, s, pressure)
450 rho = 1.0 / specvol
451 drho_dt = -dsv_dt * rho**2
452 drho_ds = -dsv_ds * rho**2
453
455
456!> Compute the in situ density of sea water (rho in [kg m-3]) and the compressibility
457!! (drho/dp = C_sound^-2, stored as drho_dp [s2 m-2]) from absolute salinity (sal [g kg-1]),
458!! conservative temperature (T [degC]), and pressure [Pa], using the specific volume
459!! polynomial fit from Roquet et al. (2015).
460elemental subroutine calculate_compress_elem_roquet_spv(this, T, S, pressure, rho, drho_dp)
461 class(roquet_spv_eos), intent(in) :: this !< This EOS
462 real, intent(in) :: t !< Conservative temperature [degC]
463 real, intent(in) :: s !< Absolute salinity [g kg-1]
464 real, intent(in) :: pressure !< pressure [Pa]
465 real, intent(out) :: rho !< In situ density [kg m-3]
466 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
467 !! (also the inverse of the square of sound speed)
468 !! [s2 m-2]
469
470 ! Local variables
471 real :: zp ! Pressure [Pa]
472 real :: zt ! Conservative temperature [degC]
473 real :: zs ! The square root of absolute salinity with an offset normalized
474 ! by an assumed salinity range [nondim]
475 real :: dsv_00p_dp ! Derivative of the pressure-dependent reference specific volume profile with
476 ! pressure [m3 kg-1 Pa-1]
477 real :: dsv_ts_dp ! Derivative of the specific volume anomaly from the reference profile with
478 ! pressure [m3 kg-1 Pa-1]
479 real :: sv_00p ! A pressure-dependent but temperature and salinity independent contribution to
480 ! specific volume at the reference temperature and salinity [m3 kg-1]
481 real :: sv_ts ! specific volume without a pressure-dependent contribution [m3 kg-1]
482 real :: sv_ts0 ! A contribution to specific volume from temperature and salinity anomalies at
483 ! the surface pressure [m3 kg-1]
484 real :: sv_ts1 ! A temperature and salinity dependent specific volume contribution that is
485 ! proportional to pressure [m3 kg-1 Pa-1]
486 real :: sv_ts2 ! A temperature and salinity dependent specific volume contribution that is
487 ! proportional to pressure**2 [m3 kg-1 Pa-2]
488 real :: sv_ts3 ! A temperature and salinity dependent specific volume contribution that is
489 ! proportional to pressure**3 [m3 kg-1 Pa-3]
490 real :: sv_0s0 ! Salinity dependent specific volume at the surface pressure and zero temperature [m3 kg-1]
491 real :: dspecvol_dp ! The partial derivative of specific volume with pressure [m3 kg-1 Pa-1]
492
493 ! The following algorithm was published by Roquet et al. (2015), intended for use
494 ! with NEMO, but it is not necessarily the algorithm used in NEMO ocean model.
495
496 ! Conversions to the units used here.
497 zt = t
498 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
499 zp = pressure
500
501 ! The next two lines should be used if it is necessary to convert potential temperature and
502 ! practical salinity to conservative temperature and absolute salinity.
503 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
504 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
505
506 sv_ts3 = spv003 + (zs*spv103 + zt*spv013)
507 sv_ts2 = spv002 + (zs*(spv102 + zs*spv202) &
508 + zt*(spv012 + (zs*spv112 + zt*spv022)) )
509 sv_ts1 = spv001 + (zs*(spv101 + zs*(spv201 + zs*(spv301 + zs*spv401))) &
510 + zt*(spv011 + (zs*(spv111 + zs*(spv211 + zs*spv311)) &
511 + zt*(spv021 + (zs*(spv121 + zs*spv221) &
512 + zt*(spv031 + (zs*spv131 + zt*spv041)) )) )) )
513
514 sv_ts0 = zt*(spv010 &
515 + (zs*(spv110 + zs*(spv210 + zs*(spv310 + zs*(spv410 + zs*spv510)))) &
516 + zt*(spv020 + (zs*(spv120 + zs*(spv220 + zs*(spv320 + zs*spv420))) &
517 + zt*(spv030 + (zs*(spv130 + zs*(spv230 + zs*spv330)) &
518 + zt*(spv040 + (zs*(spv140 + zs*spv240) &
519 + zt*(spv050 + (zs*spv150 + zt*spv060)) )) )) )) ) )
520
521 sv_0s0 = spv000 + zs*(spv100 + zs*(spv200 + zs*(spv300 + zs*(spv400 + zs*(spv500 + zs*spv600)))))
522
523 sv_00p = zp*(v00 + zp*(v01 + zp*(v02 + zp*(v03 + zp*(v04 + zp*v05)))))
524
525 sv_ts = (sv_ts0 + sv_0s0) + zp*(sv_ts1 + zp*(sv_ts2 + zp*sv_ts3))
526 ! specvol = SV_TS + SV_00p ! In situ specific volume [m3 kg-1]
527 rho = 1.0 / (sv_ts + sv_00p) ! In situ density [kg m-3]
528
529 dsv_00p_dp = v00 + zp*(2.*v01 + zp*(3.*v02 + zp*(4.*v03 + zp*(5.*v04 + zp*(6.*v05)))))
530 dsv_ts_dp = sv_ts1 + zp*(2.*sv_ts2 + zp*(3.*sv_ts3))
531 dspecvol_dp = dsv_ts_dp + dsv_00p_dp ! [m3 kg-1 Pa-1]
532 drho_dp = -dspecvol_dp * rho**2 ! Compressibility [s2 m-2]
533
535
536!> Second derivatives of specific volume with respect to temperature, salinity, and pressure for a
537!! 1-d array inputs and outputs using the specific volume polynomial fit from Roquet et al. (2015).
538elemental subroutine calc_spec_vol_second_derivs_elem_roquet_spv(T, S, P, &
539 dSV_ds_ds, dSV_ds_dt, dSV_dt_dt, dSV_ds_dp, dSV_dt_dp)
540 real, intent(in) :: t !< Conservative temperature [degC]
541 real, intent(in) :: s !< Absolute salinity [g kg-1]
542 real, intent(in) :: p !< Pressure [Pa]
543 real, intent(inout) :: dsv_ds_ds !< Second derivative of specific volume with respect
544 !! to salinity [m3 kg-1 ppt-2]
545 real, intent(inout) :: dsv_ds_dt !< Second derivative of specific volume with respect
546 !! to salinity and temperature [m3 kg-1 ppt-1 degC-1]
547 real, intent(inout) :: dsv_dt_dt !< Second derivative of specific volume with respect
548 !! to temperature [m3 kg-1 degC-2]
549 real, intent(inout) :: dsv_ds_dp !< Second derivative of specific volume with respect to pressure
550 !! and salinity [m3 kg-1 ppt-1 Pa-1]
551 real, intent(inout) :: dsv_dt_dp !< Second derivative of specific volume with respect to pressure
552 !! and temperature [m3 kg-1 degC-1 Pa-1]
553 ! Local variables
554 real :: zp ! Pressure [Pa]
555 real :: zt ! Conservative temperature [degC]
556 real :: zs ! The square root of absolute salinity with an offset normalized
557 ! by an assumed salinity range [nondim]
558 real :: i_s ! The inverse of zs [nondim]
559 real :: d2sv_p0 ! A contribution to one of the second derivatives that is independent of pressure [various]
560 real :: d2sv_p1 ! A contribution to one of the second derivatives that is proportional to pressure [various]
561 real :: d2sv_p2 ! A contribution to one of the second derivatives that is proportional to pressure**2 [various]
562 real :: d2sv_p3 ! A contribution to one of the second derivatives that is proportional to pressure**3 [various]
563
564 ! Conversions to the units used here.
565 zt = t
566 zs = sqrt( abs( s + rdeltas ) * r1_s0 ) ! square root of normalized salinity plus an offset [nondim]
567 zp = p
568
569 ! The next two lines should be used if it is necessary to convert potential temperature and
570 ! practical salinity to conservative temperature and absolute salinity.
571 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
572 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
573
574 i_s = 1.0 / zs
575
576 ! Find dSV_ds_ds
577 d2sv_p3 = -spv103*i_s**2
578 d2sv_p2 = -(spv102 + zt*spv112)*i_s**2
579 d2sv_p1 = (3.*spv301 + (zt*(3.*spv311) + zs*(8.*spv401))) &
580 - ( spv101 + zt*(spv111 + zt*(spv121 + zt*spv131)) )*i_s**2
581 d2sv_p0 = (3.*spv300 + (zs*(8.*spv400 + zs*(15.*spv500 + zs*(24.*spv600))) &
582 + zt*(3.*spv310 + (zs*(8.*spv410 + zs*(15.*spv510)) &
583 + zt*(3.*spv320 + (zs*(8.*spv420) + zt*(3.*spv330))) )) )) &
584 - (spv100 + zt*(spv110 + zt*(spv120 + zt*(spv130 + zt*(spv140 + zt*spv150)))) )*i_s**2
585 dsv_ds_ds = (0.5*r1_s0)**2 * ((d2sv_p0 + zp*(d2sv_p1 + zp*(d2sv_p2 + zp*d2sv_p3))) * i_s)
586
587 ! Find dSV_ds_dt
588 d2sv_p2 = spv112
589 d2sv_p1 = spv111 + (zs*(2.*spv211 + zs*(3.*spv311)) &
590 + zt*(2.*spv121 + (zs*(4.*spv221) + zt*(3.*spv131))) )
591 d2sv_p0 = spv110 + (zs*(2.*spv210 + zs*(3.*spv310 + zs*(4.*spv410 + zs*(5.*spv510)))) &
592 + zt*(2.*spv120 + (zs*(4.*spv220 + zs*(6.*spv320 + zs*(8.*spv420))) &
593 + zt*(3.*spv130 + (zs*(6.*spv230 + zs*(9.*spv330)) &
594 + zt*(4.*spv140 + (zs*(8.*spv240) &
595 + zt*(5.*spv150))) )) )) )
596 dsv_ds_dt = (0.5*r1_s0) * ((d2sv_p0 + zp*(d2sv_p1 + zp*d2sv_p2)) * i_s)
597
598 ! Find dSV_dt_dt
599 d2sv_p2 = 2.*spv022
600 d2sv_p1 = 2.*spv021 + (zs*(2.*spv121 + zs*(2.*spv221)) &
601 + zt*(6.*spv031 + (zs*(6.*spv131) + zt*(12.*spv041))) )
602 d2sv_p0 = 2.*spv020 + (zs*(2.*spv120 + zs*( 2.*spv220 + zs*( 2.*spv320 + zs * (2.*spv420)))) &
603 + zt*(6.*spv030 + (zs*( 6.*spv130 + zs*( 6.*spv230 + zs * (6.*spv330))) &
604 + zt*(12.*spv040 + (zs*(12.*spv140 + zs *(12.*spv240)) &
605 + zt*(20.*spv050 + (zs*(20.*spv150) &
606 + zt*(30.*spv060) )) )) )) )
607 dsv_dt_dt = d2sv_p0 + zp*(d2sv_p1 + zp*d2sv_p2)
608
609 ! Find dSV_ds_dp
610 d2sv_p2 = 3.*spv103
611 d2sv_p1 = 2.*spv102 + (zs*(4.*spv202) + zt*(2.*spv112))
612 d2sv_p0 = spv101 + (zs*(2.*spv201 + zs*(3.*spv301 + zs*(4.*spv401))) &
613 + zt*(spv111 + (zs*(2.*spv211 + zs*(3.*spv311)) &
614 + zt*( spv121 + (zs*(2.*spv221) + zt*spv131)) )) )
615 dsv_ds_dp = ((d2sv_p0 + zp*(d2sv_p1 + zp*d2sv_p2)) * i_s) * (0.5*r1_s0)
616
617 ! Find dSV_dt_dp
618 d2sv_p2 = 3.*spv013
619 d2sv_p1 = 2.*spv012 + (zs*(2.*spv112) + zt*(4.*spv022))
620 d2sv_p0 = spv011 + (zs*(spv111 + zs*( spv211 + zs* spv311)) &
621 + zt*(2.*spv021 + (zs*(2.*spv121 + zs*(2.*spv221)) &
622 + zt*(3.*spv031 + (zs*(3.*spv131) + zt*(4.*spv041))) )) )
623 dsv_dt_dp = d2sv_p0 + zp*(d2sv_p1 + zp*d2sv_p2)
624
626
627!> Second derivatives of density with respect to temperature, salinity, and pressure for a
628!! 1-d array inputs and outputs using the specific volume polynomial fit from Roquet et al. (2015).
629elemental subroutine calculate_density_second_derivs_elem_roquet_spv(this, T, S, pressure, &
630 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
631 class(roquet_spv_eos), intent(in) :: this !< This EOS
632 real, intent(in) :: t !< Conservative temperature [degC]
633 real, intent(in) :: s !< Absolute salinity [g kg-1]
634 real, intent(in) :: pressure !< Pressure [Pa]
635 real, intent(inout) :: drho_ds_ds !< Second derivative of density with respect
636 !! to salinity [kg m-3 ppt-2]
637 real, intent(inout) :: drho_ds_dt !< Second derivative of density with respect
638 !! to salinity and temperature [kg m-3 ppt-1 degC-1]
639 real, intent(inout) :: drho_dt_dt !< Second derivative of density with respect
640 !! to temperature [kg m-3 degC-2]
641 real, intent(inout) :: drho_ds_dp !< Second derivative of density with respect to pressure
642 !! and salinity [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1]
643 real, intent(inout) :: drho_dt_dp !< Second derivative of density with respect to pressure
644 !! and temperature [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
645
646 ! Local variables
647 real :: rho ! The in situ density [kg m-3]
648 real :: drho_dp ! The partial derivative of density with pressure
649 ! (also the inverse of the square of sound speed) [s2 m-2]
650 real :: dsv_dt ! The partial derivative of specific volume with
651 ! conservative temperature [m3 kg-1 degC-1]
652 real :: dsv_ds ! The partial derivative of specific volume with
653 ! absolute salinity [m3 kg-1 ppt-1]
654 real :: dsv_ds_ds ! Second derivative of specific volume with respect
655 ! to salinity [m3 kg-1 ppt-2]
656 real :: dsv_ds_dt ! Second derivative of specific volume with respect
657 ! to salinity and temperature [m3 kg-1 ppt-1 degC-1]
658 real :: dsv_dt_dt ! Second derivative of specific volume with respect
659 ! to temperature [m3 kg-1 degC-2]
660 real :: dsv_ds_dp ! Second derivative of specific volume with respect to pressure
661 ! and salinity [m3 kg-1 ppt-1 Pa-1]
662 real :: dsv_dt_dp ! Second derivative of specific volume with respect to pressure
663 ! and temperature [m3 kg-1 degC-1 Pa-1]
664
666 dsv_ds_ds, dsv_ds_dt, dsv_dt_dt, dsv_ds_dp, dsv_dt_dp)
667 call this%calculate_specvol_derivs_elem(t, s, pressure, dsv_dt, dsv_ds)
668 call this%calculate_compress_elem(t, s, pressure, rho, drho_dp)
669
670 ! Find drho_ds_ds
671 drho_ds_ds = rho**2 * (2.0*rho*dsv_ds**2 - dsv_ds_ds)
672
673 ! Find drho_ds_dt
674 drho_ds_dt = rho**2 * (2.0*rho*(dsv_dt*dsv_ds) - dsv_ds_dt)
675
676 ! Find drho_dt_dt
677 drho_dt_dt = rho**2 * (2.0*rho*dsv_dt**2 - dsv_dt_dt)
678
679 ! Find drho_ds_dp
680 drho_ds_dp = -rho * (2.0*dsv_ds * drho_dp + rho * dsv_ds_dp)
681
682 ! Find drho_dt_dp
683 drho_dt_dp = -rho * (2.0*dsv_dt * drho_dp + rho * dsv_dt_dp)
684
686
687!> Return the range of temperatures, salinities and pressures for which the Roquet et al. (2015)
688!! expression for specific volume has been fitted to observations. Care should be taken when
689!! applying this equation of state outside of its fit range.
690subroutine eos_fit_range_roquet_spv(this, T_min, T_max, S_min, S_max, p_min, p_max)
691 class(roquet_spv_eos), intent(in) :: this !< This EOS
692 real, optional, intent(out) :: T_min !< The minimum conservative temperature over which this EoS is fitted [degC]
693 real, optional, intent(out) :: T_max !< The maximum conservative temperature over which this EoS is fitted [degC]
694 real, optional, intent(out) :: S_min !< The minimum absolute salinity over which this EoS is fitted [g kg-1]
695 real, optional, intent(out) :: S_max !< The maximum absolute salinity over which this EoS is fitted [g kg-1]
696 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
697 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
698
699 if (present(t_min)) t_min = -6.0
700 if (present(t_max)) t_max = 40.0
701 if (present(s_min)) s_min = 0.0
702 if (present(s_max)) s_max = 42.0
703 if (present(p_min)) p_min = 0.0
704 if (present(p_max)) p_max = 1.0e8
705
706end subroutine eos_fit_range_roquet_spv
707
708!> Calculate the in-situ density for 1D arraya inputs and outputs.
709subroutine calculate_density_array_roquet_spv(this, T, S, pressure, rho, start, npts, rho_ref)
710 class(roquet_spv_eos), intent(in) :: this !< This EOS
711 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
712 real, dimension(:), intent(in) :: S !< Salinity [PSU]
713 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
714 real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]
715 integer, intent(in) :: start !< The starting index for calculations
716 integer, intent(in) :: npts !< The number of values to calculate
717 real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]
718
719 ! Local variables
720 integer :: j
721
722 if (present(rho_ref)) then
723 do j = start, start+npts-1
724 rho(j) = density_anomaly_elem_roquet_spv(this, t(j), s(j), pressure(j), rho_ref)
725 enddo
726 else
727 do j = start, start+npts-1
728 rho(j) = density_elem_roquet_spv(this, t(j), s(j), pressure(j))
729 enddo
730 endif
731
733
734!> Calculate the in-situ specific volume for 1D array inputs and outputs.
735subroutine calculate_spec_vol_array_roquet_spv(this, T, S, pressure, specvol, start, npts, spv_ref)
736 class(roquet_spv_eos), intent(in) :: this !< This EOS
737 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
738 real, dimension(:), intent(in) :: S !< Salinity [PSU]
739 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
740 real, dimension(:), intent(out) :: specvol !< In situ specific volume [m3 kg-1]
741 integer, intent(in) :: start !< The starting index for calculations
742 integer, intent(in) :: npts !< The number of values to calculate
743 real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
744
745 ! Local variables
746 integer :: j
747
748 if (present(spv_ref)) then
749 do j = start, start+npts-1
750 specvol(j) = spec_vol_anomaly_elem_roquet_spv(this, t(j), s(j), pressure(j), spv_ref)
751 enddo
752 else
753 do j = start, start+npts-1
754 specvol(j) = spec_vol_elem_roquet_spv(this, t(j), s(j), pressure(j) )
755 enddo
756 endif
757
759
760!> \namespace mom_eos_Roquet_SpV
761!!
762!! \section section_EOS_Roquet_SpV NEMO equation of state
763!!
764!! Fabien Roquet and colleagues developed this equation of state using a simple polynomial fit
765!! to the TEOS-10 equation of state expressions for specific, for efficiency when used with a
766!! non-Boussinesq ocean model. This particular equation of state is a balance between an
767!! accuracy that matches the TEOS-10 density to better than observational uncertainty with a
768!! polynomial form that can be evaluated quickly despite having 55 terms.
769!!
770!! \subsection section_EOS_Roquet_Spv_references References
771!!
772!! Roquet, F., Madec, G., McDougall, T. J., and Barker, P. M., 2015:
773!! Accurate polynomial expressions for the density and specific volume
774!! of seawater using the TEOS-10 standard. Ocean Modelling, 90:29-43.
775
776end module mom_eos_roquet_spv