MOM_TFreeze.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!> Freezing point expressions
6module mom_tfreeze
7
8!********+*********+*********+*********+*********+*********+*********+**
9!* The subroutines in this file determine the potential temperature *
10!* or conservative temperature at which sea-water freezes. *
11!********+*********+*********+*********+*********+*********+*********+**
12use gsw_mod_toolbox, only : gsw_ct_freezing_exact
13
14implicit none ; private
15
16public calculate_tfreeze_linear, calculate_tfreeze_millero, calculate_tfreeze_teos10
17public calculate_tfreeze_teos_poly
18
19!> Compute the freezing point potential temperature [degC] from salinity [ppt] and
20!! pressure [Pa] using a simple linear expression, with coefficients passed in as arguments.
21interface calculate_tfreeze_linear
23end interface calculate_tfreeze_linear
24
25!> Compute the freezing point potential temperature [degC] from salinity [PSU] and
26!! pressure [Pa] using the expression from Millero (1978) (and in appendix A of Gill 1982),
27!! but with the of the pressure dependence changed from 7.53e-8 to 7.75e-8 to make this an
28!! expression for potential temperature (not in situ temperature), using a
29!! value that is correct at the freezing point at 35 PSU and 5e6 Pa (500 dbar).
30interface calculate_tfreeze_millero
32end interface calculate_tfreeze_millero
33
34!> Compute the freezing point conservative temperature [degC] from absolute salinity [g kg-1]
35!! and pressure [Pa] using the TEOS10 package.
36interface calculate_tfreeze_teos10
38end interface calculate_tfreeze_teos10
39
40!> Compute the freezing point conservative temperature [degC] from absolute salinity [g kg-1] and
41!! pressure [Pa] using a rescaled and refactored version of the expressions from the TEOS10 package.
42interface calculate_tfreeze_teos_poly
44end interface calculate_tfreeze_teos_poly
45
46contains
47
48!> This subroutine computes the freezing point potential temperature [degC] from
49!! salinity [ppt], and pressure [Pa] using a simple linear expression,
50!! with coefficients passed in as arguments.
51subroutine calculate_tfreeze_linear_scalar(S, pres, T_Fr, TFr_S0_P0, &
52 dTFr_dS, dTFr_dp)
53 real, intent(in) :: S !< salinity [ppt].
54 real, intent(in) :: pres !< pressure [Pa].
55 real, intent(out) :: T_Fr !< Freezing point potential temperature [degC].
56 real, intent(in) :: TFr_S0_P0 !< The freezing point at S=0, p=0 [degC].
57 real, intent(in) :: dTFr_dS !< The derivative of freezing point with salinity,
58 !! [degC ppt-1].
59 real, intent(in) :: dTFr_dp !< The derivative of freezing point with pressure,
60 !! [degC Pa-1].
61
62 t_fr = (tfr_s0_p0 + dtfr_ds*s) + dtfr_dp*pres
63
65
66!> This subroutine computes an array of freezing point potential temperatures
67!! [degC] from salinity [ppt], and pressure [Pa] using a simple
68!! linear expression, with coefficients passed in as arguments.
69subroutine calculate_tfreeze_linear_array(S, pres, T_Fr, start, npts, &
70 TFr_S0_P0, dTFr_dS, dTFr_dp)
71 real, dimension(:), intent(in) :: S !< salinity [ppt].
72 real, dimension(:), intent(in) :: pres !< pressure [Pa].
73 real, dimension(:), intent(out) :: T_Fr !< Freezing point potential temperature [degC].
74 integer, intent(in) :: start !< the starting point in the arrays.
75 integer, intent(in) :: npts !< the number of values to calculate.
76 real, intent(in) :: TFr_S0_P0 !< The freezing point at S=0, p=0, [degC].
77 real, intent(in) :: dTFr_dS !< The derivative of freezing point with salinity,
78 !! [degC ppt-1].
79 real, intent(in) :: dTFr_dp !< The derivative of freezing point with pressure,
80 !! [degC Pa-1].
81 integer :: j
82
83 do j=start,start+npts-1
84 t_fr(j) = (tfr_s0_p0 + dtfr_ds*s(j)) + dtfr_dp*pres(j)
85 enddo
86
88
89!> This subroutine computes the freezing point potential temperature
90!! [degC] from salinity [ppt], and pressure [Pa] using the expression
91!! from Millero (1978) (and in appendix A of Gill 1982), but with the of the
92!! pressure dependence changed from 7.53e-8 to 7.75e-8 to make this an
93!! expression for potential temperature (not in situ temperature), using a
94!! value that is correct at the freezing point at 35 PSU and 5e6 Pa (500 dbar).
95subroutine calculate_tfreeze_millero_scalar(S, pres, T_Fr)
96 real, intent(in) :: S !< Salinity [PSU]
97 real, intent(in) :: pres !< Pressure [Pa]
98 real, intent(out) :: T_Fr !< Freezing point potential temperature [degC]
99
100 ! Local variables
101 real, parameter :: cS1 = -0.0575 ! A term in the freezing point fit [degC PSU-1]
102 real, parameter :: cS3_2 = 1.710523e-3 ! A term in the freezing point fit [degC PSU-3/2]
103 real, parameter :: cS2 = -2.154996e-4 ! A term in the freezing point fit [degC PSU-2]
104 real, parameter :: dTFr_dp = -7.75e-8 ! Derivative of freezing point with pressure [degC Pa-1]
105
106 t_fr = s*(cs1 + (cs3_2 * sqrt(max(s, 0.0)) + cs2 * s)) + dtfr_dp*pres
107
109
110!> This subroutine computes the freezing point potential temperature
111!! [degC] from salinity [ppt], and pressure [Pa] using the expression
112!! from Millero (1978) (and in appendix A of Gill 1982), but with the
113!! pressure dependence changed from 7.53e-8 to 7.75e-8 to make this an
114!! expression for potential temperature (not in situ temperature), using a
115!! value that is correct at the freezing point at 35 PSU and 5e6 Pa (500 dbar).
116subroutine calculate_tfreeze_millero_array(S, pres, T_Fr, start, npts)
117 real, dimension(:), intent(in) :: S !< Salinity [PSU].
118 real, dimension(:), intent(in) :: pres !< Pressure [Pa].
119 real, dimension(:), intent(out) :: T_Fr !< Freezing point potential temperature [degC].
120 integer, intent(in) :: start !< The starting point in the arrays.
121 integer, intent(in) :: npts !< The number of values to calculate.
122
123 ! Local variables
124 real, parameter :: cS1 = -0.0575 ! A term in the freezing point fit [degC PSU-1]
125 real, parameter :: cS3_2 = 1.710523e-3 ! A term in the freezing point fit [degC PSU-3/2]
126 real, parameter :: cS2 = -2.154996e-4 ! A term in the freezing point fit [degC PSU-2]
127 real, parameter :: dTFr_dp = -7.75e-8 ! Derivative of freezing point with pressure [degC Pa-1]
128 integer :: j
129
130 do j=start,start+npts-1
131 t_fr(j) = s(j)*(cs1 + (cs3_2 * sqrt(max(s(j), 0.0)) + cs2 * s(j))) + &
132 dtfr_dp*pres(j)
133 enddo
134
136
137!> This subroutine computes the freezing point conservative temperature [degC]
138!! from absolute salinity [g kg-1], and pressure [Pa] using a rescaled and
139!! refactored version of the polynomial expressions from the TEOS10 package.
140subroutine calculate_tfreeze_teos_poly_scalar(S, pres, T_Fr)
141 real, intent(in) :: S !< Absolute salinity [g kg-1].
142 real, intent(in) :: pres !< Pressure [Pa].
143 real, intent(out) :: T_Fr !< Freezing point conservative temperature [degC].
144
145 ! Local variables
146 real, dimension(1) :: S0 ! Salinity at a point [g kg-1]
147 real, dimension(1) :: pres0 ! Pressure at a point [Pa]
148 real, dimension(1) :: tfr0 ! The freezing temperature [degC]
149
150 s0(1) = s
151 pres0(1) = pres
152
153 call calculate_tfreeze_teos_poly_array(s0, pres0, tfr0, 1, 1)
154 t_fr = tfr0(1)
155
157
158!> This subroutine computes the freezing point conservative temperature [degC]
159!! from absolute salinity [g kg-1], and pressure [Pa] using a rescaled and
160!! refactored version of the polynomial expressions from the TEOS10 package.
161subroutine calculate_tfreeze_teos_poly_array(S, pres, T_Fr, start, npts)
162 real, dimension(:), intent(in) :: S !< absolute salinity [g kg-1].
163 real, dimension(:), intent(in) :: pres !< Pressure [Pa].
164 real, dimension(:), intent(out) :: T_Fr !< Freezing point conservative temperature [degC].
165 integer, intent(in) :: start !< The starting point in the arrays
166 integer, intent(in) :: npts !< The number of values to calculate
167
168 ! Local variables
169 real :: Sa ! Absolute salinity [g kg-1] = [ppt]
170 real :: rS ! Square root of salinity [ppt1/2]
171 ! The coefficients here use the notation TFab for contributions proportional to S**a/2 * P**b.
172 real, parameter :: TF00 = 0.017947064327968736 ! Freezing point coefficient [degC]
173 real, parameter :: TF20 = -6.076099099929818e-2 ! Freezing point coefficient [degC ppt-1]
174 real, parameter :: TF30 = 4.883198653547851e-3 ! Freezing point coefficient [degC ppt-3/2]
175 real, parameter :: TF40 = -1.188081601230542e-3 ! Freezing point coefficient [degC ppt-2]
176 real, parameter :: TF50 = 1.334658511480257e-4 ! Freezing point coefficient [degC ppt-5/2]
177 real, parameter :: TF60 = -8.722761043208607e-6 ! Freezing point coefficient [degC ppt-3]
178 real, parameter :: TF70 = 2.082038908808201e-7 ! Freezing point coefficient [degC ppt-7/2]
179 real, parameter :: TF01 = -7.389420998107497e-8 ! Freezing point coefficient [degC Pa-1]
180 real, parameter :: TF21 = -9.891538123307282e-11 ! Freezing point coefficient [degC ppt-1 Pa-1]
181 real, parameter :: TF31 = -8.987150128406496e-13 ! Freezing point coefficient [degC ppt-3/2 Pa-1]
182 real, parameter :: TF41 = 1.054318231187074e-12 ! Freezing point coefficient [degC ppt-2 Pa-1]
183 real, parameter :: TF51 = 3.850133554097069e-14 ! Freezing point coefficient [degC ppt-5/2 Pa-1]
184 real, parameter :: TF61 = -2.079022768390933e-14 ! Freezing point coefficient [degC ppt-3 Pa-1]
185 real, parameter :: TF71 = 1.242891021876471e-15 ! Freezing point coefficient [degC ppt-7/2 Pa-1]
186 real, parameter :: TF02 = -2.110913185058476e-16 ! Freezing point coefficient [degC Pa-2]
187 real, parameter :: TF22 = 3.831132432071728e-19 ! Freezing point coefficient [degC ppt-1 Pa-2]
188 real, parameter :: TF32 = 1.065556599652796e-19 ! Freezing point coefficient [degC ppt-3/2 Pa-2]
189 real, parameter :: TF42 = -2.078616693017569e-20 ! Freezing point coefficient [degC ppt-2 Pa-2]
190 real, parameter :: TF52 = 1.596435439942262e-21 ! Freezing point coefficient [degC ppt-5/2 Pa-2]
191 real, parameter :: TF03 = 2.295491578006229e-25 ! Freezing point coefficient [degC Pa-3]
192 real, parameter :: TF23 = -7.997496801694032e-27 ! Freezing point coefficient [degC ppt-1 Pa-3]
193 real, parameter :: TF33 = 8.756340772729538e-28 ! Freezing point coefficient [degC ppt-3/2 Pa-3]
194 real, parameter :: TF43 = 1.338002171109174e-29 ! Freezing point coefficient [degC ppt-2 Pa-3]
195 integer :: j
196
197 do j=start,start+npts-1
198 rs = sqrt(max(s(j), 0.0))
199 t_fr(j) = (tf00 + s(j)*(tf20 + rs*(tf30 + rs*(tf40 + rs*(tf50 + rs*(tf60 + rs*tf70)))))) &
200 + pres(j)*( (tf01 + s(j)*(tf21 + rs*(tf31 + rs*(tf41 + rs*(tf51 + rs*(tf61 + rs*tf71)))))) &
201 + pres(j)*((tf02 + s(j)*(tf22 + rs*(tf32 + rs*(tf42 + rs* tf52)))) &
202 + pres(j)*(tf03 + s(j)*(tf23 + rs*(tf33 + rs* tf43))) ) )
203 enddo
204
206
207!> This subroutine computes the freezing point conservative temperature [degC]
208!! from absolute salinity [g kg-1], and pressure [Pa] using the
209!! TEOS10 package.
210subroutine calculate_tfreeze_teos10_scalar(S, pres, T_Fr)
211 real, intent(in) :: S !< Absolute salinity [g kg-1].
212 real, intent(in) :: pres !< Pressure [Pa].
213 real, intent(out) :: T_Fr !< Freezing point conservative temperature [degC].
214
215 ! Local variables
216 real, dimension(1) :: S0 ! Salinity at a point [g kg-1]
217 real, dimension(1) :: pres0 ! Pressure at a point [Pa]
218 real, dimension(1) :: tfr0 ! The freezing temperature [degC]
219
220 s0(1) = s
221 pres0(1) = pres
222
223 call calculate_tfreeze_teos10_array(s0, pres0, tfr0, 1, 1)
224 t_fr = tfr0(1)
225
227
228!> This subroutine computes the freezing point conservative temperature [degC]
229!! from absolute salinity [g kg-1], and pressure [Pa] using the
230!! TEOS10 package.
231subroutine calculate_tfreeze_teos10_array(S, pres, T_Fr, start, npts)
232 real, dimension(:), intent(in) :: S !< absolute salinity [g kg-1].
233 real, dimension(:), intent(in) :: pres !< pressure [Pa].
234 real, dimension(:), intent(out) :: T_Fr !< Freezing point conservative temperature [degC].
235 integer, intent(in) :: start !< the starting point in the arrays.
236 integer, intent(in) :: npts !< the number of values to calculate.
237
238 ! Local variables
239 real, parameter :: Pa2db = 1.e-4 ! The conversion factor from Pa to dbar [dbar Pa-1]
240 real :: zp ! Pressures in [dbar]
241 integer :: j
242 ! Assume sea-water contains no dissolved air.
243 real, parameter :: saturation_fraction = 0.0 ! Air saturation fraction in seawater [nondim]
244
245 do j=start,start+npts-1
246 !Conversions
247 zp = pres(j)* pa2db !Convert pressure from Pascal to decibar
248
249 if (s(j) < -1.0e-10) cycle !Can we assume safely that this is a missing value?
250 t_fr(j) = gsw_ct_freezing_exact(s(j), zp, saturation_fraction)
251 enddo
252
254
255end module mom_tfreeze