← back to index

src/equation_of_state/MOM_EOS_Roquet_rho.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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
6module MOM_EOS_Roquet_rho
7
8use MOM_EOS_base_type, only : EOS_base
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_density_array_2d for efficiency
185 procedure :: calculate_density_array_2d => calculate_density_array_2d_Roquet_rho
186 !> Local implementation of generic calculate_density_array_3d for efficiency
187 procedure :: calculate_density_array_3d => calculate_density_array_3d_Roquet_rho
188 !> Local implementation of generic calculate_spec_vol_array for efficiency
189 procedure :: calculate_spec_vol_array => calculate_spec_vol_array_Roquet_rho
190 !> Local implementation of generic calculate_density_derivs_2d for efficiency
191 procedure :: calculate_density_derivs_2d => calculate_density_derivs_2d_Roquet_rho
192 !> Local implementation of generic calculate_density_derivs_3d for efficiency
193 procedure :: calculate_density_derivs_3d => calculate_density_derivs_3d_Roquet_rho
194 !> Local implementation of generic calculate_density_second_derivs_2d for efficiency
195 procedure :: calculate_density_second_derivs_2d => calculate_density_second_derivs_2d_Roquet_rho
196
197end type Roquet_rho_EOS
198
199contains
200
201!> In situ density of sea water from Roquet et al., 2015 [kg m-3]
202!!
203!! This is an elemental function that can be applied to any combination of scalar and array inputs.
20412119572real elemental function density_elem_Roquet_rho_loc(T, S, pressure)
205 real, intent(in) :: T !< Conservative temperature [degC]
206 real, intent(in) :: S !< Absolute salinity [g kg-1]
207 real, intent(in) :: pressure !< Pressure [Pa]
208
209 ! Local variables
210 real :: zp ! Pressure [Pa]
211 real :: zt ! Conservative temperature [degC]
212 real :: zs ! The square root of absolute salinity with an offset normalized
213 ! by an assumed salinity range [nondim]
214 real :: rho00p ! A pressure-dependent but temperature and salinity independent contribution to
215 ! density at the reference temperature and salinity [kg m-3]
216 real :: rhoTS ! Density without a pressure-dependent contribution [kg m-3]
217 real :: rhoTS0 ! A contribution to density from temperature and salinity anomalies at the
218 ! surface pressure [kg m-3]
219 real :: rhoTS1 ! A density contribution proportional to pressure [kg m-3 Pa-1]
220 real :: rhoTS2 ! A density contribution proportional to pressure**2 [kg m-3 Pa-2]
221 real :: rhoTS3 ! A density contribution proportional to pressure**3 [kg m-3 Pa-3]
222 real :: rho0S0 ! Salinity dependent density at the surface pressure and zero temperature [kg m-3]
223
224 ! The following algorithm was published by Roquet et al. (2015), intended for use with NEMO.
225
226 ! Conversions to the units used here.
22712119572 zt = T
22812119572 zs = SQRT( ABS( S + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim]
22912119572 zp = pressure
230
231 ! The next two lines should be used if it is necessary to convert potential temperature and
232 ! practical salinity to conservative temperature and absolute salinity.
233 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
234 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
235
23612119572 rhoTS3 = EOS003 + (zs*EOS103 + zt*EOS013)
237 rhoTS2 = EOS002 + (zs*(EOS102 + zs*EOS202) &
23812119572 + zt*(EOS012 + (zs*EOS112 + zt*EOS022)) )
239 rhoTS1 = EOS001 + (zs*(EOS101 + zs*(EOS201 + zs*(EOS301 + zs*EOS401))) &
240 + zt*(EOS011 + (zs*(EOS111 + zs*(EOS211 + zs*EOS311)) &
241 + zt*(EOS021 + (zs*(EOS121 + zs*EOS221) &
24212119572 + zt*(EOS031 + (zs*EOS131 + zt*EOS041)) )) )) )
243 rhoTS0 = zt*(EOS010 &
244 + (zs*(EOS110 + zs*(EOS210 + zs*(EOS310 + zs*(EOS410 + zs*EOS510)))) &
245 + zt*(EOS020 + (zs*(EOS120 + zs*(EOS220 + zs*(EOS320 + zs*EOS420))) &
246 + zt*(EOS030 + (zs*(EOS130 + zs*(EOS230 + zs*EOS330)) &
247 + zt*(EOS040 + (zs*(EOS140 + zs*EOS240) &
24812119572 + zt*(EOS050 + (zs*EOS150 + zt*EOS060)) )) )) )) ) )
249
25012119572 rho0S0 = EOS000 + zs*(EOS100 + zs*(EOS200 + zs*(EOS300 + zs*(EOS400 + zs*(EOS500 + zs*EOS600)))))
251
25212119572 rho00p = zp*(R00 + zp*(R01 + zp*(R02 + zp*(R03 + zp*(R04 + zp*R05)))))
253
25412119572 rhoTS = (rhoTS0 + rho0S0) + zp*(rhoTS1 + zp*(rhoTS2 + zp*rhoTS3))
25512119572 density_elem_Roquet_rho_loc = rhoTS + rho00p ! In situ density [kg m-3]
256
25712119572end function density_elem_Roquet_rho_loc
258
259!> Wrapper for density_elem_Roquet_rho_loc created to preserve API while calling
260!! density_elem_Roquet_rho without "this" variable that causes runtime errors on
261!! gpu runs with nvfortran.
2626480000real elemental function density_elem_Roquet_rho(this, T, S, pressure)
263 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
264 real, intent(in) :: T !< Conservative temperature [degC]
265 real, intent(in) :: S !< Absolute salinity [g kg-1]
266 real, intent(in) :: pressure !< Pressure [Pa]
267
2686480000 density_elem_Roquet_rho = density_elem_Roquet_rho_loc(T, S, pressure)
2696480000end function density_elem_Roquet_rho
270
271!> In situ density anomaly of sea water from Roquet et al., 2015 [kg m-3]
272!!
273!! This is an elemental function that can be applied to any combination of scalar and array inputs.
274461736000real elemental function density_anomaly_elem_Roquet_rho_loc(T, S, pressure, rho_ref)
275 real, intent(in) :: T !< Conservative temperature [degC]
276 real, intent(in) :: S !< Absolute salinity [g kg-1]
277 real, intent(in) :: pressure !< Pressure [Pa]
278 real, intent(in) :: rho_ref !< A reference density [kg m-3]
279
280 ! Local variables
281 real :: zp ! Pressure [Pa]
282 real :: zt ! Conservative temperature [degC]
283 real :: zs ! The square root of absolute salinity with an offset normalized
284 ! by an assumed salinity range [nondim]
285 real :: rho00p ! A pressure-dependent but temperature and salinity independent contribution to
286 ! density at the reference temperature and salinity [kg m-3]
287 real :: rhoTS ! Density without a pressure-dependent contribution [kg m-3]
288 real :: rhoTS0 ! A contribution to density from temperature and salinity anomalies at the
289 ! surface pressure [kg m-3]
290 real :: rhoTS1 ! A density contribution proportional to pressure [kg m-3 Pa-1]
291 real :: rhoTS2 ! A density contribution proportional to pressure**2 [kg m-3 Pa-2]
292 real :: rhoTS3 ! A density contribution proportional to pressure**3 [kg m-3 Pa-3]
293 real :: rho0S0 ! Salinity dependent density at the surface pressure and zero temperature [kg m-3]
294
295 ! The following algorithm was published by Roquet et al. (2015), intended for use with NEMO.
296
297 ! Conversions to the units used here.
298461736000 zt = T
299461736000 zs = SQRT( ABS( S + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim]
300461736000 zp = pressure
301
302 ! The next two lines should be used if it is necessary to convert potential temperature and
303 ! practical salinity to conservative temperature and absolute salinity.
304 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
305 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
306
307461736000 rhoTS3 = EOS003 + (zs*EOS103 + zt*EOS013)
308 rhoTS2 = EOS002 + (zs*(EOS102 + zs*EOS202) &
309461736000 + zt*(EOS012 + (zs*EOS112 + zt*EOS022)) )
310 rhoTS1 = EOS001 + (zs*(EOS101 + zs*(EOS201 + zs*(EOS301 + zs*EOS401))) &
311 + zt*(EOS011 + (zs*(EOS111 + zs*(EOS211 + zs*EOS311)) &
312 + zt*(EOS021 + (zs*(EOS121 + zs*EOS221) &
313461736000 + zt*(EOS031 + (zs*EOS131 + zt*EOS041)) )) )) )
314 rhoTS0 = zt*(EOS010 &
315 + (zs*(EOS110 + zs*(EOS210 + zs*(EOS310 + zs*(EOS410 + zs*EOS510)))) &
316 + zt*(EOS020 + (zs*(EOS120 + zs*(EOS220 + zs*(EOS320 + zs*EOS420))) &
317 + zt*(EOS030 + (zs*(EOS130 + zs*(EOS230 + zs*EOS330)) &
318 + zt*(EOS040 + (zs*(EOS140 + zs*EOS240) &
319461736000 + zt*(EOS050 + (zs*EOS150 + zt*EOS060)) )) )) )) ) )
320
321461736000 rho0S0 = EOS000 + zs*(EOS100 + zs*(EOS200 + zs*(EOS300 + zs*(EOS400 + zs*(EOS500 + zs*EOS600)))))
322
323461736000 rho00p = zp*(R00 + zp*(R01 + zp*(R02 + zp*(R03 + zp*(R04 + zp*R05)))))
324
325461736000 rho0S0 = rho0S0 - rho_ref
326
327461736000 rhoTS = (rhoTS0 + rho0S0) + zp*(rhoTS1 + zp*(rhoTS2 + zp*rhoTS3))
328461736000 density_anomaly_elem_Roquet_rho_loc = rhoTS + rho00p ! In situ density [kg m-3]
329
330461736000end function density_anomaly_elem_Roquet_rho_loc
331
332!> Wrapper for density_anomaly_elem_Roquet_rho_loc created to preserve API while calling
333!! density_anomaly_elem_Roquet_rho_loc without "this" variable that causes runtime errors on
334!! gpu runs with nvfortran.
3350real elemental function density_anomaly_elem_Roquet_rho(this, T, S, pressure, rho_ref)
336 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
337 real, intent(in) :: T !< Conservative temperature [degC]
338 real, intent(in) :: S !< Absolute salinity [g kg-1]
339 real, intent(in) :: pressure !< Pressure [Pa]
340 real, intent(in) :: rho_ref !< A reference density [kg m-3]
341
3420 density_anomaly_elem_Roquet_rho = density_anomaly_elem_Roquet_rho_loc(T, S, pressure, rho_ref)
3430end function density_anomaly_elem_Roquet_rho
344
345!> In situ specific volume of sea water from Roquet et al., 2015 [kg m-3]
346!!
347!! This is an elemental function that can be applied to any combination of scalar and array inputs.
3480real elemental function spec_vol_elem_Roquet_rho(this, T, S, pressure)
349 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
350 real, intent(in) :: T !< Conservative temperature [degC]
351 real, intent(in) :: S !< Absolute salinity [g kg-1]
352 real, intent(in) :: pressure !< Pressure [Pa]
353
3540 spec_vol_elem_Roquet_rho = 1. / density_elem_Roquet_rho(this, T, S, pressure)
355
3560end function spec_vol_elem_Roquet_rho
357
358!> In situ specific volume anomaly of sea water from Roquet et al., 2015 [kg m-3]
359!!
360!! This is an elemental function that can be applied to any combination of scalar and array inputs.
3610real elemental function spec_vol_anomaly_elem_Roquet_rho(this, T, S, pressure, spv_ref)
362 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
363 real, intent(in) :: T !< Conservative temperature [degC]
364 real, intent(in) :: S !< Absolute salinity [g kg-1]
365 real, intent(in) :: pressure !< Pressure [Pa]
366 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
367
3680 spec_vol_anomaly_elem_Roquet_rho = 1. / density_elem_Roquet_rho(this, T, S, pressure)
3690 spec_vol_anomaly_elem_Roquet_rho = spec_vol_anomaly_elem_Roquet_rho - spv_ref
370
3710end function spec_vol_anomaly_elem_Roquet_rho
372
373!> For a given thermodynamic state, calculate the derivatives of density with conservative
374!! temperature and absolute salinity, using the density polynomial fit EOS from Roquet et al. (2015).
37570323852elemental subroutine calculate_density_derivs_elem_Roquet_rho_loc(T, S, pressure, drho_dT, drho_dS)
376 real, intent(in) :: T !< Conservative temperature [degC]
377 real, intent(in) :: S !< Absolute salinity [g kg-1]
378 real, intent(in) :: pressure !< Pressure [Pa]
379 real, intent(out) :: drho_dT !< The partial derivative of density with potential
380 !! temperature [kg m-3 degC-1]
381 real, intent(out) :: drho_dS !< The partial derivative of density with salinity,
382 !! in [kg m-3 ppt-1]
383
384 ! Local variables
385 real :: zp ! Pressure [Pa]
386 real :: zt ! Conservative temperature [degC]
387 real :: zs ! The square root of absolute salinity with an offset normalized
388 ! by an assumed salinity range [nondim]
389 real :: dRdzt0 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1]
390 ! from temperature anomalies at the surface pressure
391 real :: dRdzt1 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1 Pa-1]
392 ! proportional to pressure
393 real :: dRdzt2 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1 Pa-2]
394 ! proportional to pressure**2
395 real :: dRdzt3 ! A contribution to the partial derivative of density with temperature [kg m-3 degC-1 Pa-3]
396 ! proportional to pressure**3
397 real :: dRdzs0 ! A contribution to the partial derivative of density with
398 ! salinity [kg m-3 ppt-1] from temperature anomalies at the surface pressure
399 real :: dRdzs1 ! A contribution to the partial derivative of density with
400 ! salinity [kg m-3 ppt-1 Pa-1] proportional to pressure
401 real :: dRdzs2 ! A contribution to the partial derivative of density with
402 ! salinity [kg m-3 ppt-1 Pa-2] proportional to pressure**2
403 real :: dRdzs3 ! A contribution to the partial derivative of density with
404 ! salinity [kg m-3 ppt-1 Pa-3] proportional to pressure**3
405
406 ! Conversions to the units used here.
40770323852 zt = T
40870323852 zs = SQRT( ABS( S + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim]
40970323852 zp = pressure
410
411 ! The next two lines should be used if it is necessary to convert potential temperature and
412 ! practical salinity to conservative temperature and absolute salinity.
413 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
414 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
415
416 ! Find the partial derivative of density with temperature
41770323852 dRdzt3 = ALP003
41870323852 dRdzt2 = ALP002 + (zs*ALP102 + zt*ALP012)
419 dRdzt1 = ALP001 + (zs*(ALP101 + zs*(ALP201 + zs*ALP301)) &
420 + zt*(ALP011 + (zs*(ALP111 + zs*ALP211) &
42170323852 + zt*(ALP021 + (zs*ALP121 + zt*ALP031)) )) )
422 dRdzt0 = ALP000 + (zs*(ALP100 + zs*(ALP200 + zs*(ALP300 + zs*(ALP400 + zs*ALP500)))) &
423 + zt*(ALP010 + (zs*(ALP110 + zs*(ALP210 + zs*(ALP310 + zs*ALP410))) &
424 + zt*(ALP020 + (zs*(ALP120 + zs*(ALP220 + zs*ALP320)) &
425 + zt*(ALP030 + (zt*(ALP040 + (zs*ALP140 + zt*ALP050)) &
42670323852 + zs*(ALP130 + zs*ALP230) )) )) )) )
427
42870323852 drho_dT = dRdzt0 + zp*(dRdzt1 + zp*(dRdzt2 + zp*dRdzt3))
429
430 ! Find the partial derivative of density with salinity
43170323852 dRdzs3 = BET003
43270323852 dRdzs2 = BET002 + (zs*BET102 + zt*BET012)
433 dRdzs1 = BET001 + (zs*(BET101 + zs*(BET201 + zs*BET301)) &
434 + zt*(BET011 + (zs*(BET111 + zs*BET211) &
43570323852 + zt*(BET021 + (zs*BET121 + zt*BET031)) )) )
436 dRdzs0 = BET000 + (zs*(BET100 + zs*(BET200 + zs*(BET300 + zs*(BET400 + zs*BET500)))) &
437 + zt*(BET010 + (zs*(BET110 + zs*(BET210 + zs*(BET310 + zs*BET410))) &
438 + zt*(BET020 + (zs*(BET120 + zs*(BET220 + zs*BET320)) &
439 + zt*(BET030 + (zt*(BET040 + (zs*BET140 + zt*BET050)) &
44070323852 + zs*(BET130 + zs*BET230) )) )) )) )
441
442 ! The division by zs here is because zs = sqrt(S + S0), so drho_dS = dzs_dS * drho_dzs = (0.5 / zs) * drho_dzs
44370323852 drho_dS = (dRdzs0 + zp*(dRdzs1 + zp*(dRdzs2 + zp * dRdzs3))) / zs
444
44570323852end subroutine calculate_density_derivs_elem_Roquet_rho_loc
446
447!> Wrapper for calculate_density_derivs_elem_Roquet_rho_loc created to preserve API while
448!! calling calculate_density_derivs_elem_Roquet_rho without "this" variable that causes
449!! runtime errors on gpu runs with nvfortran.
45043118172elemental subroutine calculate_density_derivs_elem_Roquet_rho(this, T, S, pressure, drho_dT, drho_dS)
451 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
452 real, intent(in) :: T !< Conservative temperature [degC]
453 real, intent(in) :: S !< Absolute salinity [g kg-1]
454 real, intent(in) :: pressure !< Pressure [Pa]
455 real, intent(out) :: drho_dT !< The partial derivative of density with potential
456 !! temperature [kg m-3 degC-1]
457 real, intent(out) :: drho_dS !< The partial derivative of density with salinity,
458 !! in [kg m-3 ppt-1]
459
46043118172 call calculate_density_derivs_elem_Roquet_rho_loc(T, S, pressure, drho_dT, drho_dS)
461
46243118172end subroutine calculate_density_derivs_elem_Roquet_rho
463
464!> Second derivatives of density with respect to temperature, salinity, and pressure.
465!! Free-function form without "this" for use in do concurrent GPU regions with nvfortran.
4660elemental subroutine calculate_density_second_derivs_elem_Roquet_rho_loc(T, S, pressure, &
467 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
468 real, intent(in) :: T !< Conservative temperature [degC]
469 real, intent(in) :: S !< Absolute salinity [g kg-1]
470 real, intent(in) :: pressure !< Pressure [Pa]
471 real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
472 !! to S [kg m-3 ppt-2]
473 real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
474 !! to T [kg m-3 ppt-1 degC-1]
475 real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
476 !! to T [kg m-3 degC-2]
477 real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
478 !! to pressure [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1]
479 real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
480 !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
481
482 ! Local variables
483 real :: zp ! Pressure [Pa]
484 real :: zt ! Conservative temperature [degC]
485 real :: zs ! The square root of absolute salinity with an offset normalized
486 ! by an assumed salinity range [nondim]
487 real :: I_s ! The inverse of zs [nondim]
488 real :: d2R_p0 ! A contribution to one of the second derivatives that is independent of pressure [various]
489 real :: d2R_p1 ! A contribution to one of the second derivatives that is proportional to pressure [various]
490 real :: d2R_p2 ! A contribution to one of the second derivatives that is proportional to pressure**2 [various]
491 real :: d2R_p3 ! A contribution to one of the second derivatives that is proportional to pressure**3 [various]
492
493 ! Conversions to the units used here.
4940 zt = T
4950 zs = SQRT( ABS( S + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim]
4960 zp = pressure
497
498 ! The next two lines should be used if it is necessary to convert potential temperature and
499 ! practical salinity to conservative temperature and absolute salinity.
500 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
501 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
502
5030 I_s = 1.0 / zs
504
505 ! Find drho_ds_ds
5060 d2R_p3 = -EOS103*I_s**2
5070 d2R_p2 = -(EOS102 + zt*EOS112)*I_s**2
508 d2R_p1 = (3.*EOS301 + (zt*(3.*EOS311) + zs*(8.*EOS401))) &
5090 - ( EOS101 + zt*(EOS111 + zt*(EOS121 + zt*EOS131)) )*I_s**2
510 d2R_p0 = (3.*EOS300 + (zs*(8.*EOS400 + zs*(15.*EOS500 + zs*(24.*EOS600))) &
511 + zt*(3.*EOS310 + (zs*(8.*EOS410 + zs*(15.*EOS510)) &
512 + zt*(3.*EOS320 + (zs*(8.*EOS420) + zt*(3.*EOS330))) )) )) &
5130 - (EOS100 + zt*(EOS110 + zt*(EOS120 + zt*(EOS130 + zt*(EOS140 + zt*EOS150)))) )*I_s**2
5140 drho_dS_dS = (0.5*r1_S0)**2 * ((d2R_p0 + zp*(d2R_p1 + zp*(d2R_p2 + zp*d2R_p3))) * I_s)
515
516 ! Find drho_ds_dt
5170 d2R_p2 = EOS112
518 d2R_p1 = EOS111 + (zs*(2.*EOS211 + zs*(3.*EOS311)) &
5190 + zt*(2.*EOS121 + (zs*(4.*EOS221) + zt*(3.*EOS131))) )
520 d2R_p0 = EOS110 + (zs*(2.*EOS210 + zs*(3.*EOS310 + zs*(4.*EOS410 + zs*(5.*EOS510)))) &
521 + zt*(2.*EOS120 + (zs*(4.*EOS220 + zs*(6.*EOS320 + zs*(8.*EOS420))) &
522 + zt*(3.*EOS130 + (zs*(6.*EOS230 + zs*(9.*EOS330)) &
523 + zt*(4.*EOS140 + (zs*(8.*EOS240) &
5240 + zt*(5.*EOS150))) )) )) )
5250 drho_ds_dt = (0.5*r1_S0) * ((d2R_p0 + zp*(d2R_p1 + zp*d2R_p2)) * I_s)
526
527 ! Find drho_dt_dt
5280 d2R_p2 = 2.*EOS022
529 d2R_p1 = 2.*EOS021 + (zs*(2.*EOS121 + zs*(2.*EOS221)) &
5300 + zt*(6.*EOS031 + (zs*(6.*EOS131) + zt*(12.*EOS041))) )
531 d2R_p0 = 2.*EOS020 + (zs*(2.*EOS120 + zs*( 2.*EOS220 + zs*( 2.*EOS320 + zs * (2.*EOS420)))) &
532 + zt*(6.*EOS030 + (zs*( 6.*EOS130 + zs*( 6.*EOS230 + zs * (6.*EOS330))) &
533 + zt*(12.*EOS040 + (zs*(12.*EOS140 + zs *(12.*EOS240)) &
534 + zt*(20.*EOS050 + (zs*(20.*EOS150) &
5350 + zt*(30.*EOS060) )) )) )) )
5360 drho_dt_dt = (d2R_p0 + zp*(d2R_p1 + zp*d2R_p2))
537
538 ! Find drho_ds_dp
5390 d2R_p2 = 3.*EOS103
5400 d2R_p1 = 2.*EOS102 + (zs*(4.*EOS202) + zt*(2.*EOS112))
541 d2R_p0 = EOS101 + (zs*(2.*EOS201 + zs*(3.*EOS301 + zs*(4.*EOS401))) &
542 + zt*(EOS111 + (zs*(2.*EOS211 + zs*(3.*EOS311)) &
5430 + zt*( EOS121 + (zs*(2.*EOS221) + zt*EOS131)) )) )
5440 drho_ds_dp = ((d2R_p0 + zp*(d2R_p1 + zp*d2R_p2)) * I_s) * (0.5*r1_S0)
545
546 ! Find drho_dt_dp
5470 d2R_p2 = 3.*EOS013
5480 d2R_p1 = 2.*EOS012 + (zs*(2.*EOS112) + zt*(4.*EOS022))
549 d2R_p0 = EOS011 + (zs*(EOS111 + zs*( EOS211 + zs* EOS311)) &
550 + zt*(2.*EOS021 + (zs*(2.*EOS121 + zs*(2.*EOS221)) &
5510 + zt*(3.*EOS031 + (zs*(3.*EOS131) + zt*(4.*EOS041))) )) )
5520 drho_dt_dp = (d2R_p0 + zp*(d2R_p1 + zp*d2R_p2))
553
5540end subroutine calculate_density_second_derivs_elem_Roquet_rho_loc
555
556!> Wrapper for calculate_density_second_derivs_elem_Roquet_rho_loc created to preserve API
557!! while calling without "this" variable that causes runtime errors on GPU with nvfortran.
5580elemental subroutine calculate_density_second_derivs_elem_Roquet_rho(this, T, S, pressure, &
559 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
560 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
561 real, intent(in) :: T !< Conservative temperature [degC]
562 real, intent(in) :: S !< Absolute salinity [g kg-1]
563 real, intent(in) :: pressure !< Pressure [Pa]
564 real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
565 !! to S [kg m-3 ppt-2]
566 real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
567 !! to T [kg m-3 ppt-1 degC-1]
568 real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
569 !! to T [kg m-3 degC-2]
570 real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
571 !! to pressure [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1]
572 real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
573 !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
574
575 call calculate_density_second_derivs_elem_Roquet_rho_loc(T, S, pressure, &
5760 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
577
5780end subroutine calculate_density_second_derivs_elem_Roquet_rho
579
580!> Calculate the partial derivatives of specific volume with temperature and salinity
581!! using the density polynomial fit EOS from Roquet et al. (2015).
5826480000elemental subroutine calculate_specvol_derivs_elem_Roquet_rho(this, T, S, pressure, dSV_dT, dSV_dS)
583 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
584 real, intent(in) :: T !< Conservative temperature [degC]
585 real, intent(in) :: S !< Absolute salinity [g kg-1]
586 real, intent(in) :: pressure !< Pressure [Pa]
587 real, intent(inout) :: dSV_dT !< The partial derivative of specific volume with
588 !! potential temperature [m3 kg-1 degC-1]
589 real, intent(inout) :: dSV_dS !< The partial derivative of specific volume with
590 !! salinity [m3 kg-1 ppt-1]
591 ! Local variables
592 real :: rho ! In situ density [kg m-3]
593 real :: dRho_dT ! Derivative of density with temperature [kg m-3 degC-1]
594 real :: dRho_dS ! Derivative of density with salinity [kg m-3 ppt-1]
595
5966480000 call this%calculate_density_derivs_elem(T, S, pressure, drho_dT, drho_dS)
5976480000 rho = this%density_elem(T, S, pressure)
5986480000 dSV_dT = -dRho_DT/(rho**2)
5996480000 dSV_dS = -dRho_DS/(rho**2)
600
6016480000end subroutine calculate_specvol_derivs_elem_Roquet_rho
602
603!> Compute the in situ density of sea water (rho in [kg m-3]) and the compressibility
604!! (drho/dp = C_sound^-2, stored as drho_dp [s2 m-2]) from absolute salinity (sal [g kg-1]),
605!! conservative temperature (T [degC]), and pressure [Pa], using the density polynomial
606!! fit EOS from Roquet et al. (2015).
6070elemental subroutine calculate_compress_elem_Roquet_rho(this, T, S, pressure, rho, drho_dp)
608 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
609 real, intent(in) :: T !< Conservative temperature [degC]
610 real, intent(in) :: S !< Absolute salinity [g kg-1]
611 real, intent(in) :: pressure !< Pressure [Pa]
612 real, intent(out) :: rho !< In situ density [kg m-3]
613 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
614 !! (also the inverse of the square of sound speed)
615 !! [s2 m-2]
616 ! Local variables
617 real :: zp ! Pressure [Pa]
618 real :: zt ! Conservative temperature [degC]
619 real :: zs ! The square root of absolute salinity with an offset normalized
620 ! by an assumed salinity range [nondim]
621 real :: drho00p_dp ! Derivative of the pressure-dependent reference density profile with pressure [kg m-3 Pa-1]
622 real :: drhoTS_dp ! Derivative of the density anomaly from the reference profile with pressure [kg m-3 Pa-1]
623 real :: rho00p ! The pressure-dependent (but temperature and salinity independent) reference
624 ! density profile [kg m-3]
625 real :: rhoTS ! Density anomaly from the reference profile [kg m-3]
626 real :: rhoTS0 ! A contribution to density from temperature and salinity anomalies at the
627 ! surface pressure [kg m-3]
628 real :: rhoTS1 ! A density contribution proportional to pressure [kg m-3 Pa-1]
629 real :: rhoTS2 ! A density contribution proportional to pressure**2 [kg m-3 Pa-2]
630 real :: rhoTS3 ! A density contribution proportional to pressure**3 [kg m-3 Pa-3]
631 real :: rho0S0 ! Salinity dependent density at the surface pressure and zero temperature [kg m-3]
632
633 ! The following algorithm was published by Roquet et al. (2015), intended for use with NEMO.
634 ! Conversions to the units used here.
6350 zt = T
6360 zs = SQRT( ABS( S + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim]
6370 zp = pressure
638
639 ! The next two lines should be used if it is necessary to convert potential temperature and
640 ! practical salinity to conservative temperature and absolute salinity.
641 ! zt = gsw_ct_from_pt(S,T) ! Convert potential temp to conservative temp [degC]
642 ! zs = SQRT( ABS( gsw_sr_from_sp(S) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity.
643
6440 rhoTS3 = EOS003 + (zs*EOS103 + zt*EOS013)
645 rhoTS2 = EOS002 + (zs*(EOS102 + zs*EOS202) &
6460 + zt*(EOS012 + (zs*EOS112 + zt*EOS022)) )
647 rhoTS1 = EOS001 + (zs*(EOS101 + zs*(EOS201 + zs*(EOS301 + zs*EOS401))) &
648 + zt*(EOS011 + (zs*(EOS111 + zs*(EOS211 + zs*EOS311)) &
649 + zt*(EOS021 + (zs*(EOS121 + zs*EOS221) &
6500 + zt*(EOS031 + (zs*EOS131 + zt*EOS041)) )) )) )
651
652 rhoTS0 = zt*(EOS010 &
653 + (zs*(EOS110 + zs*(EOS210 + zs*(EOS310 + zs*(EOS410 + zs*EOS510)))) &
654 + zt*(EOS020 + (zs*(EOS120 + zs*(EOS220 + zs*(EOS320 + zs*EOS420))) &
655 + zt*(EOS030 + (zs*(EOS130 + zs*(EOS230 + zs*EOS330)) &
656 + zt*(EOS040 + (zs*(EOS140 + zs*EOS240) &
6570 + zt*(EOS050 + (zs*EOS150 + zt*EOS060)) )) )) )) ) )
658
6590 rho0S0 = EOS000 + zs*(EOS100 + zs*(EOS200 + zs*(EOS300 + zs*(EOS400 + zs*(EOS500 + zs*EOS600)))))
660
6610 rho00p = zp*(R00 + zp*(R01 + zp*(R02 + zp*(R03 + zp*(R04 + zp*R05)))))
662
6630 rhoTS = (rhoTS0 + rho0S0) + zp*(rhoTS1 + zp*(rhoTS2 + zp*rhoTS3))
6640 rho = rhoTS + rho00p ! In situ density [kg m-3]
665
6660 drho00p_dp = R00 + zp*(2.*R01 + zp*(3.*R02 + zp*(4.*R03 + zp*(5.*R04 + zp*(6.*R05)))))
6670 drhoTS_dp = rhoTS1 + zp*(2.*rhoTS2 + zp*(3.*rhoTS3))
6680 drho_dp = drhoTS_dp + drho00p_dp ! Compressibility [s2 m-2]
669
6700end subroutine calculate_compress_elem_Roquet_rho
671
672!> Return the range of temperatures, salinities and pressures for which the Roquet et al. (2015)
673!! expression for in situ density has been fitted to observations. Care should be taken when
674!! applying this equation of state outside of its fit range.
6750subroutine EoS_fit_range_Roquet_rho(this, T_min, T_max, S_min, S_max, p_min, p_max)
676 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
677 real, optional, intent(out) :: T_min !< The minimum conservative temperature over which this EoS is fitted [degC]
678 real, optional, intent(out) :: T_max !< The maximum conservative temperature over which this EoS is fitted [degC]
679 real, optional, intent(out) :: S_min !< The minimum absolute salinity over which this EoS is fitted [g kg-1]
680 real, optional, intent(out) :: S_max !< The maximum absolute salinity over which this EoS is fitted [g kg-1]
681 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
682 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
683
6840 if (present(T_min)) T_min = -6.0
6850 if (present(T_max)) T_max = 40.0
6860 if (present(S_min)) S_min = 0.0
6870 if (present(S_max)) S_max = 42.0
6880 if (present(p_min)) p_min = 0.0
6890 if (present(p_max)) p_max = 1.0e8
690
6910end subroutine EoS_fit_range_Roquet_rho
692
693!> Calculate the in-situ density for 1D arraya inputs and outputs.
69444738subroutine calculate_density_array_Roquet_rho(this, T, S, pressure, rho, start, npts, rho_ref)
695 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
696 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
697 real, dimension(:), intent(in) :: S !< Salinity [PSU]
698 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
699 real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]
700 integer, intent(in) :: start !< The starting index for calculations
701 integer, intent(in) :: npts !< The number of values to calculate
702 real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]
703
704 ! Local variables
705 integer :: j
706
70744738 if (present(rho_ref)) then
7080 do j = start, start+npts-1
7090 rho(j) = density_anomaly_elem_Roquet_rho(this, T(j), S(j), pressure(j), rho_ref)
710 enddo
711 else
7125502774 do j = start, start+npts-1
7135502774 rho(j) = density_elem_Roquet_rho_loc(T(j), S(j), pressure(j))
714 enddo
715 endif
716
71744738end subroutine calculate_density_array_Roquet_rho
718
719!> Calculate the in-situ density for 2D array inputs and outputs.
72024subroutine calculate_density_array_2d_Roquet_rho(this, T, S, pressure, rho, dom, rho_ref)
721 class(Roquet_rho_EOS), intent(in) :: this
722 !< This EOS
723 real, intent(in) :: T(:,:)
724 !< Conservative temperature [degC]
725 real, intent(in) :: S(:,:)
726 !< Absolute salinity [g kg-1]
727 real, intent(in) :: pressure(:,:)
728 !< Pressure [Pa]
729 real, intent(out) :: rho(:,:)
730 !< In situ density [kg m-3]
731 integer, intent(in) :: dom(2,2)
732 !< Index bounds of domain. First index is rank, second is bounds
733 real, optional, intent(in) :: rho_ref
734 !< A reference density [kg m-3]
735
736 integer :: is, ie, js, je
737 integer :: i, j
738
73924 is = dom(1,1) ; ie = dom(1,2)
74024 js = dom(2,1) ; je = dom(2,2)
741
74224 if (present(rho_ref)) then
7430 do concurrent (j=js:je, i=is:ie)
744 rho(i,j) = density_anomaly_elem_Roquet_rho_loc(T(i,j), S(i,j), &
7450 pressure(i,j), rho_ref)
746 enddo
747 else
74824 do concurrent (j=js:je, i=is:ie)
749184488 rho(i,j) = density_elem_Roquet_rho_loc( T(i,j), S(i,j), pressure(i,j))
750 enddo
751 endif
75224end subroutine calculate_density_array_2d_Roquet_rho
753
754!> Calculate the in-situ density for 3D array inputs and outputs.
755329400subroutine calculate_density_array_3d_Roquet_rho(this, T, S, pressure, rho, dom, rho_ref)
756 class(Roquet_rho_EOS), intent(in) :: this
757 !< This EOS
758 real, intent(in) :: T(:,:,:)
759 !< Conservative temperature [degC]
760 real, intent(in) :: S(:,:,:)
761 !< Absolute salinity [g kg-1]
762 real, intent(in) :: pressure(:,:,:)
763 !< Pressure [Pa]
764 real, intent(out) :: rho(:,:,:)
765 !< In situ density [kg m-3]
766 integer, intent(in) :: dom(3,2)
767 !< Index bounds of domain. First index is rank, second is bounds
768 real, optional, intent(in) :: rho_ref
769 !< A reference density [kg m-3]
770
771 integer :: is, ie, js, je, ks, ke
772 integer :: i, j, k
773
774329400 is = dom(1,1) ; ie = dom(1,2)
775329400 js = dom(2,1) ; je = dom(2,2)
776329400 ks = dom(3,1) ; ke = dom(3,2)
777
778 ! The element functions are called via their free-function (_loc) forms rather than
779 ! through the polymorphic "this" binding, which causes runtime errors in do concurrent
780 ! regions offloaded to the GPU with nvfortran.
781329400 if (present(rho_ref)) then
782329400 do concurrent (k=ks:ke, j=js:je, i=is:ie)
783 rho(i,j,k) = density_anomaly_elem_Roquet_rho_loc(T(i,j,k), S(i,j,k), &
7841385537400 pressure(i,j,k), rho_ref)
785 enddo
786 else
7870 do concurrent (k=ks:ke, j=js:je, i=is:ie)
7880 rho(i,j,k) = density_elem_Roquet_rho_loc(T(i,j,k), S(i,j,k), pressure(i,j,k))
789 enddo
790 endif
791329400end subroutine calculate_density_array_3d_Roquet_rho
792
793!> Calculate the in-situ density derivatives for 2D array inputs and outputs.
79424subroutine calculate_density_derivs_2d_Roquet_rho(this, T, S, pressure, &
79524 drho_dT, drho_dS, dom)
796 class(Roquet_rho_EOS), intent(in) :: this
797 !< This EOS
798 real, intent(in) :: T(:,:)
799 !< Conservative temperature [degC]
800 real, intent(in) :: S(:,:)
801 !< Absolute salinity [g kg-1]
802 real, intent(in) :: pressure(:,:)
803 !< Pressure [Pa]
804 real, intent(out) :: drho_dT(:,:)
805 !< Partial derivative of density with potential temperature [kg m-3 degC-1]
806 real, intent(out) :: drho_dS(:,:)
807 !< Partial derivative of density with salinity [kg m-3 ppt-1]
808 integer, intent(in) :: dom(2,2)
809 !< Index bounds of domain. First index is rank, second is bounds
810
811 integer :: is, ie, js, je
812 integer :: i, j
813
81424 is = dom(1,1) ; ie = dom(1,2)
81524 js = dom(2,1) ; je = dom(2,2)
816
817 ! NOTE: There is an implicit copy of `this` which cannot yet be prevented.
818
81924 do concurrent (j=js:je, i=is:ie)
820 call calculate_density_derivs_elem_Roquet_rho_loc(T(i,j), S(i,j), &
821177876 pressure(i,j), drho_dT(i,j), drho_dS(i,j))
822 enddo
82324end subroutine calculate_density_derivs_2d_Roquet_rho
824
825!> Calculate the in-situ density derivatives for 3D array inputs and outputs.
826223776subroutine calculate_density_derivs_3d_Roquet_rho(this, T, S, pressure, &
827223776 drho_dT, drho_dS, dom)
828 class(Roquet_rho_EOS), intent(in) :: this
829 !< This EOS
830 real, intent(in) :: T(:,:,:)
831 !< Conservative temperature [degC]
832 real, intent(in) :: S(:,:,:)
833 !< Absolute salinity [g kg-1]
834 real, intent(in) :: pressure(:,:,:)
835 !< Pressure [Pa]
836 real, intent(out) :: drho_dT(:,:,:)
837 !< Partial derivative of density with potential temperature [kg m-3 degC-1]
838 real, intent(out) :: drho_dS(:,:,:)
839 !< Partial derivative of density with salinity [kg m-3 ppt-1]
840 integer, intent(in) :: dom(3,2)
841 !< Index bounds of domain. First index is rank, second is bounds
842
843 integer :: is, ie, js, je, ks, ke
844 integer :: i, j, k
845
846223776 is = dom(1,1) ; ie = dom(1,2)
847223776 js = dom(2,1) ; je = dom(2,2)
848223776 ks = dom(3,1) ; ke = dom(3,2)
849
850 ! The element subroutine is called via its free-function (_loc) form rather than
851 ! through the polymorphic "this" binding, which causes runtime errors in do concurrent
852 ! regions offloaded to the GPU with nvfortran.
853223776 do concurrent (k=ks:ke, j=js:je, i=is:ie)
854 call calculate_density_derivs_elem_Roquet_rho_loc(T(i,j,k), S(i,j,k), &
85568098944 pressure(i,j,k), drho_dT(i,j,k), drho_dS(i,j,k))
856 enddo
857223776end subroutine calculate_density_derivs_3d_Roquet_rho
858
859!> Calculate the second derivatives of density for 2D array inputs and outputs.
8600subroutine calculate_density_second_derivs_2d_Roquet_rho(this, T, S, pressure, &
8610 drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, dom)
862 class(Roquet_rho_EOS), intent(in) :: this
863 !< This EOS
864 real, intent(in) :: T(:,:)
865 !< Conservative temperature [degC]
866 real, intent(in) :: S(:,:)
867 !< Absolute salinity [g kg-1]
868 real, intent(in) :: pressure(:,:)
869 !< Pressure [Pa]
870 real, intent(inout) :: drho_dS_dS(:,:)
871 !< Partial derivative of beta with respect to S [kg m-3 ppt-2]
872 real, intent(inout) :: drho_dS_dT(:,:)
873 !< Partial derivative of beta with respect to T [kg m-3 ppt-1 degC-1]
874 real, intent(inout) :: drho_dT_dT(:,:)
875 !< Partial derivative of alpha with respect to T [kg m-3 degC-2]
876 real, intent(inout) :: drho_dS_dP(:,:)
877 !< Partial derivative of beta with respect to pressure [kg m-3 ppt-1 Pa-1]
878 real, intent(inout) :: drho_dT_dP(:,:)
879 !< Partial derivative of alpha with respect to pressure [kg m-3 degC-1 Pa-1]
880 integer, intent(in) :: dom(2,2)
881 !< Index bounds of domain. First index is rank, second is bounds
882
883 integer :: is, ie, js, je
884 integer :: i, j
885
8860 is = dom(1,1) ; ie = dom(1,2)
8870 js = dom(2,1) ; je = dom(2,2)
888
889 ! The element subroutine is called via its free-function (_loc) form rather than
890 ! through the polymorphic "this" binding, which causes runtime errors in do concurrent
891 ! regions offloaded to the GPU with nvfortran.
8920 do concurrent (j=js:je, i=is:ie)
893 call calculate_density_second_derivs_elem_Roquet_rho_loc(T(i,j), S(i,j), pressure(i,j), &
8940 drho_dS_dS(i,j), drho_dS_dT(i,j), drho_dT_dT(i,j), drho_dS_dP(i,j), drho_dT_dP(i,j))
895 enddo
8960end subroutine calculate_density_second_derivs_2d_Roquet_rho
897
898!> Calculate the in-situ specific volume for 1D array inputs and outputs.
8990subroutine calculate_spec_vol_array_Roquet_rho(this, T, S, pressure, specvol, start, npts, spv_ref)
900 class(Roquet_rho_EOS), intent(in) :: this !< This EOS
901 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
902 real, dimension(:), intent(in) :: S !< Salinity [PSU]
903 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
904 real, dimension(:), intent(out) :: specvol !< In situ specific volume [m3 kg-1]
905 integer, intent(in) :: start !< The starting index for calculations
906 integer, intent(in) :: npts !< The number of values to calculate
907 real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
908
909 ! Local variables
910 integer :: j
911
9120 if (present(spv_ref)) then
9130 do j = start, start+npts-1
9140 specvol(j) = spec_vol_anomaly_elem_Roquet_rho(this, T(j), S(j), pressure(j), spv_ref)
915 enddo
916 else
9170 do j = start, start+npts-1
9180 specvol(j) = spec_vol_elem_Roquet_rho(this, T(j), S(j), pressure(j) )
919 enddo
920 endif
921
9220end subroutine calculate_spec_vol_array_Roquet_rho
923
924!> \namespace mom_eos_Roquet_rho
925!!
926!! \section section_EOS_Roquet_rho Roquet_rho equation of state
927!!
928!! Fabien Roquet and colleagues developed this equation of state using a simple polynomial fit
929!! to the TEOS-10 equation of state, for efficiency when used in the NEMO ocean model. Fabien
930!! Roquet also graciously provided the MOM6 team with the original code implementing this
931!! equation of state, although it has since been modified and extended to have capabilities
932!! mirroring those available with other equations of state in MOM6. This particular equation
933!! of state is a balance between an accuracy that matches the TEOS-10 density to better than
934!! observational uncertainty with a polynomial form that can be evaluated quickly despite having
935!! 52 terms.
936!!
937!! \subsection section_EOS_Roquet_rho_references References
938!!
939!! Roquet, F., Madec, G., McDougall, T. J., and Barker, P. M., 2015:
940!! Accurate polynomial expressions for the density and specific volume
941!! of seawater using the TEOS-10 standard. Ocean Modelling, 90:29-43.
942
9430end module MOM_EOS_Roquet_rho