MOM_EOS_TEOS10.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> The equation of state using the TEOS10 expressions
7
8use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp, gsw_ct_from_pt
9use gsw_mod_toolbox, only : gsw_rho, gsw_specvol
10use gsw_mod_toolbox, only : gsw_rho_first_derivatives, gsw_specvol_first_derivatives
11use gsw_mod_toolbox, only : gsw_rho_second_derivatives
13
14implicit none ; private
15
16public gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp, gsw_ct_from_pt
17public teos10_eos
18
19real, parameter :: pa2db = 1.e-4 !< The conversion factor from Pa to dbar [dbar Pa-1]
20
21!> The EOS_base implementation of the TEOS10 equation of state
22type, extends (eos_base) :: teos10_eos
23
24contains
25 !> Implementation of the in-situ density as an elemental function [kg m-3]
26 procedure :: density_elem => density_elem_teos10
27 !> Implementation of the in-situ density anomaly as an elemental function [kg m-3]
28 procedure :: density_anomaly_elem => density_anomaly_elem_teos10
29 !> Implementation of the in-situ specific volume as an elemental function [m3 kg-1]
30 procedure :: spec_vol_elem => spec_vol_elem_teos10
31 !> Implementation of the in-situ specific volume anomaly as an elemental function [m3 kg-1]
32 procedure :: spec_vol_anomaly_elem => spec_vol_anomaly_elem_teos10
33 !> Implementation of the calculation of derivatives of density
34 procedure :: calculate_density_derivs_elem => calculate_density_derivs_elem_teos10
35 !> Implementation of the calculation of second derivatives of density
36 procedure :: calculate_density_second_derivs_elem => calculate_density_second_derivs_elem_teos10
37 !> Implementation of the calculation of derivatives of specific volume
38 procedure :: calculate_specvol_derivs_elem => calculate_specvol_derivs_elem_teos10
39 !> Implementation of the calculation of compressibility
40 procedure :: calculate_compress_elem => calculate_compress_elem_teos10
41 !> Implementation of the range query function
42 procedure :: eos_fit_range => eos_fit_range_teos10
43
44end type teos10_eos
45
46contains
47
48!> GSW in situ density [kg m-3]
49real elemental function density_elem_teos10(this, t, s, pressure)
50 class(teos10_eos), intent(in) :: this !< This EOS
51 real, intent(in) :: t !< Conservative temperature [degC].
52 real, intent(in) :: s !< Absolute salinity [g kg-1].
53 real, intent(in) :: pressure !< pressure [Pa].
54
55 density_elem_teos10 = gsw_rho(s, t, pressure * pa2db)
56
57end function density_elem_teos10
58
59!> GSW in situ density anomaly [kg m-3]
60real elemental function density_anomaly_elem_teos10(this, t, s, pressure, rho_ref)
61 class(teos10_eos), intent(in) :: this !< This EOS
62 real, intent(in) :: t !< Conservative temperature [degC].
63 real, intent(in) :: s !< Absolute salinity [g kg-1].
64 real, intent(in) :: pressure !< pressure [Pa].
65 real, intent(in) :: rho_ref !< A reference density [kg m-3].
66
67 density_anomaly_elem_teos10 = gsw_rho(s, t, pressure * pa2db)
69
71
72!> GSW in situ specific volume [m3 kg-1]
73real elemental function spec_vol_elem_teos10(this, t, s, pressure)
74 class(teos10_eos), intent(in) :: this !< This EOS
75 real, intent(in) :: t !< Conservative temperature [degC].
76 real, intent(in) :: s !< Absolute salinity [g kg-1].
77 real, intent(in) :: pressure !< pressure [Pa].
78
79 spec_vol_elem_teos10 = gsw_specvol(s, t, pressure * pa2db)
80
81end function spec_vol_elem_teos10
82
83!> GSW in situ specific volume anomaly [m3 kg-1]
84real elemental function spec_vol_anomaly_elem_teos10(this, t, s, pressure, spv_ref)
85 class(teos10_eos), intent(in) :: this !< This EOS
86 real, intent(in) :: t !< Conservative temperature [degC].
87 real, intent(in) :: s !< Absolute salinity [g kg-1].
88 real, intent(in) :: pressure !< pressure [Pa].
89 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
90
91 spec_vol_anomaly_elem_teos10 = gsw_specvol(s, t, pressure * pa2db) - spv_ref
92
94
95!> For a given thermodynamic state, calculate the derivatives of density with conservative
96!! temperature and absolute salinity, using the TEOS10 expressions.
97elemental subroutine calculate_density_derivs_elem_teos10(this, T, S, pressure, drho_dT, drho_dS)
98 class(teos10_eos), intent(in) :: this !< This EOS
99 real, intent(in) :: t !< Conservative temperature [degC]
100 real, intent(in) :: s !< Absolute salinity [g kg-1] = [ppt]
101 real, intent(in) :: pressure !< Pressure [Pa]
102 real, intent(out) :: drho_dt !< The partial derivative of density with conservative
103 !! temperature [kg m-3 degC-1]
104 real, intent(out) :: drho_ds !< The partial derivative of density with salinity,
105 !! in [kg m-3 ppt-1]
106 ! Local variables
107 real :: zs ! Absolute salinity [g kg-1]
108 real :: zt ! Conservative temperature [degC]
109 real :: zp ! Pressure converted to decibars [dbar]
110
111 ! Conversions
112 zs = s
113 zt = t
114 zp = pressure * pa2db ! Convert pressure from Pascal to decibar
115 ! The following conversions are unnecessary because the arguments are already the right variables.
116 ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity
117 ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp
118
119 call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_ds, drho_dct=drho_dt)
120
122
123!> Calculate the 5 second derivatives of the equation of state for scalar inputs
124elemental subroutine calculate_density_second_derivs_elem_teos10(this, T, S, pressure, &
125 drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP)
126 class(teos10_eos), intent(in) :: this !< This EOS
127 real, intent(in) :: t !< Conservative temperature [degC]
128 real, intent(in) :: s !< Absolute salinity [g kg-1] = [ppt]
129 real, intent(in) :: pressure !< Pressure [Pa]
130 real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
131 !! to S [kg m-3 ppt-2]
132 real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect
133 !! to T [kg m-3 ppt-1 degC-1]
134 real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
135 !! to T [kg m-3 degC-2]
136 real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
137 !! to pressure [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1]
138 real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
139 !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1]
140 ! Local variables
141 real :: zs ! Absolute salinity [g kg-1]
142 real :: zt ! Conservative temperature [degC]
143 real :: zp ! Pressure converted to decibars [dbar]
144
145 ! Conversions
146 zs = s
147 zt = t
148 zp = pressure * pa2db ! Convert pressure from Pascal to decibar
149 ! The following conversions are unnecessary because the arguments are already the right variables.
150 ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity
151 ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp
152
153 call gsw_rho_second_derivatives(zs, zt, zp, rho_sa_sa=drho_ds_ds, rho_sa_ct=drho_ds_dt, &
154 rho_ct_ct=drho_dt_dt, rho_sa_p=drho_ds_dp, rho_ct_p=drho_dt_dp)
155
157
158!> For a given thermodynamic state, calculate the derivatives of specific volume with conservative
159!! temperature and absolute salinity, using the TEOS10 expressions.
160elemental subroutine calculate_specvol_derivs_elem_teos10(this, T, S, pressure, dSV_dT, dSV_dS)
161 class(teos10_eos), intent(in) :: this !< This EOS
162 real, intent(in) :: t !< Conservative temperature [degC]
163 real, intent(in) :: s !< Absolute salinity [g kg-1] = [ppt]
164 real, intent(in) :: pressure !< Pressure [Pa]
165 real, intent(inout) :: dsv_dt !< The partial derivative of specific volume with
166 !! conservative temperature [m3 kg-1 degC-1]
167 real, intent(inout) :: dsv_ds !< The partial derivative of specific volume with
168 !! absolute salinity [m3 kg-1 ppt-1]
169 ! Local variables
170 real :: zs ! Absolute salinity [g kg-1]
171 real :: zt ! Conservative temperature [degC]
172 real :: zp ! Pressure converted to decibars [dbar]
173
174 ! Conversions
175 zs = s
176 zt = t
177 zp = pressure * pa2db ! Convert pressure from Pascal to decibar
178 ! The following conversions are unnecessary because the arguments are already the right variables.
179 ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity
180 ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp
181
182 call gsw_specvol_first_derivatives(zs, zt, zp, v_sa=dsv_ds, v_ct=dsv_dt)
183
185
186!> This subroutine computes the in situ density of sea water (rho in
187!! [kg m-3]) and the compressibility (drho/dp = C_sound^-2)
188!! (drho_dp [s2 m-2]) from absolute salinity (sal [g kg-1]),
189!! conservative temperature (T [degC]), and pressure [Pa]. It uses the
190!! subroutines from TEOS10 website
191elemental subroutine calculate_compress_elem_teos10(this, T, S, pressure, rho, drho_dp)
192 class(teos10_eos), intent(in) :: this !< This EOS
193 real, intent(in) :: t !< Conservative temperature [degC]
194 real, intent(in) :: s !< Absolute salinity [g kg-1]
195 real, intent(in) :: pressure !< Pressure [Pa]
196 real, intent(out) :: rho !< In situ density [kg m-3]
197 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
198 !! (also the inverse of the square of sound speed)
199 !! [s2 m-2]
200
201 ! Local variables
202 real :: zs ! Absolute salinity [g kg-1]
203 real :: zt ! Conservative temperature [degC]
204 real :: zp ! Pressure converted to decibars [dbar]
205
206 ! Conversions
207 zs = s
208 zt = t
209 zp = pressure * pa2db ! Convert pressure from Pascal to decibar
210 ! The following conversions are unnecessary because the arguments are already the right variables.
211 ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity
212 ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp
213
214 rho = gsw_rho(zs, zt, zp)
215 call gsw_rho_first_derivatives(zs, zt, zp, drho_dp=drho_dp)
216
218
219!> Return the range of temperatures, salinities and pressures for which the TEOS-10
220!! equation of state has been fitted to observations. Care should be taken when
221!! applying this equation of state outside of its fit range.
222subroutine eos_fit_range_teos10(this, T_min, T_max, S_min, S_max, p_min, p_max)
223 class(teos10_eos), intent(in) :: this !< This EOS
224 real, optional, intent(out) :: T_min !< The minimum conservative temperature over which this EoS is fitted [degC]
225 real, optional, intent(out) :: T_max !< The maximum conservative temperature over which this EoS is fitted [degC]
226 real, optional, intent(out) :: S_min !< The minimum absolute salinity over which this EoS is fitted [g kg-1]
227 real, optional, intent(out) :: S_max !< The maximum absolute salinity over which this EoS is fitted [g kg-1]
228 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
229 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
230
231 if (present(t_min)) t_min = -6.0
232 if (present(t_max)) t_max = 40.0
233 if (present(s_min)) s_min = 0.0
234 if (present(s_max)) s_max = 42.0
235 if (present(p_min)) p_min = 0.0
236 if (present(p_max)) p_max = 1.0e8
237
238end subroutine eos_fit_range_teos10
239
240!> \namespace mom_eos_teos10
241!!
242!! \section section_EOS_TEOS10 TEOS10 equation of state
243!!
244!! The TEOS10 equation of state is implemented via the GSW toolbox. We recommend using the
245!! Roquet et al. forms of this equation of state.
246
247end module mom_eos_teos10