MOM_temperature_convert.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!> Functions to convert between conservative and potential temperature
7
8implicit none ; private
9
11
12!>@{ Parameters in the temperature conversion code
13real, parameter :: sprac_sref = (35.0/35.16504) ! The TEOS 10 conversion factor to go from
14 ! reference salinity to practical salinity [nondim]
15real, parameter :: i_s0 = 0.025*sprac_sref ! The inverse of a plausible range of oceanic salinities [kg g-1]
16real, parameter :: i_ts = 0.025 ! The inverse of a plausible range of oceanic temperatures [degC-1]
17real, parameter :: i_cp0 = 1.0/3991.86795711963 ! The inverse of the "specific heat" for use
18 ! with Conservative Temperature, as defined with TEOS10 [degC kg J-1]
19
20! The following are coefficients of contributions to conservative temperature as a function of the square root
21! of normalized absolute salinity with an offset (zS) and potential temperature (T) with a contribution
22! Hab * zS**a * T**b. The numbers here are copied directly from the corresponding gsw module, but
23! the expressions here do not use the same nondimensionalization for pressure or temperature as they do.
24
25real, parameter :: h00 = 61.01362420681071*i_cp0 ! Tp to Tc fit constant [degC]
26real, parameter :: h01 = 168776.46138048015*(i_cp0*i_ts) ! Tp to Tc fit T coef. [nondim]
27real, parameter :: h02 = -2735.2785605119625*(i_cp0*i_ts**2) ! Tp to Tc fit T**2 coef. [degC-1]
28real, parameter :: h03 = 2574.2164453821433*(i_cp0*i_ts**3) ! Tp to Tc fit T**3 coef. [degC-2]
29real, parameter :: h04 = -1536.6644434977543*(i_cp0*i_ts**4) ! Tp to Tc fit T**4 coef. [degC-3]
30real, parameter :: h05 = 545.7340497931629*(i_cp0*i_ts**5) ! Tp to Tc fit T**5 coef. [degC-4]
31real, parameter :: h06 = -50.91091728474331*(i_cp0*i_ts**6) ! Tp to Tc fit T**6 coef. [degC-5]
32real, parameter :: h07 = -18.30489878927802*(i_cp0*i_ts**7) ! Tp to Tc fit T**7 coef. [degC-6]
33real, parameter :: h20 = 268.5520265845071*i_cp0 ! Tp to Tc fit zS**2 coef. [degC]
34real, parameter :: h21 = -12019.028203559312*(i_cp0*i_ts) ! Tp to Tc fit zS**2 * T coef. [nondim]
35real, parameter :: h22 = 3734.858026725145*(i_cp0*i_ts**2) ! Tp to Tc fit zS**2 * T**2 coef. [degC-1]
36real, parameter :: h23 = -2046.7671145057618*(i_cp0*i_ts**3) ! Tp to Tc fit zS**2 * T**3 coef. [degC-2]
37real, parameter :: h24 = 465.28655623826234*(i_cp0*i_ts**4) ! Tp to Tc fit zS**2 * T**4 coef. [degC-3]
38real, parameter :: h25 = -0.6370820302376359*(i_cp0*i_ts**5) ! Tp to Tc fit zS**2 * T**5 coef. [degC-4]
39real, parameter :: h26 = -10.650848542359153*(i_cp0*i_ts**6) ! Tp to Tc fit zS**2 * T**6 coef. [degC-5]
40real, parameter :: h30 = 937.2099110620707*i_cp0 ! Tp to Tc fit zS**3 coef. [degC]
41real, parameter :: h31 = 588.1802812170108*(i_cp0*i_ts) ! Tp to Tc fit zS** 3* T coef. [nondim]
42real, parameter :: h32 = 248.39476522971285*(i_cp0*i_ts**2) ! Tp to Tc fit zS**3 * T**2 coef. [degC-1]
43real, parameter :: h33 = -3.871557904936333*(i_cp0*i_ts**3) ! Tp to Tc fit zS**3 * T**3 coef. [degC-2]
44real, parameter :: h34 = -2.6268019854268356*(i_cp0*i_ts**4) ! Tp to Tc fit zS**3 * T**4 coef. [degC-3]
45real, parameter :: h40 = -1687.914374187449*i_cp0 ! Tp to Tc fit zS**4 coef. [degC]
46real, parameter :: h41 = 936.3206544460336*(i_cp0*i_ts) ! Tp to Tc fit zS**4 * T coef. [nondim]
47real, parameter :: h42 = -942.7827304544439*(i_cp0*i_ts**2) ! Tp to Tc fit zS**4 * T**2 coef. [degC-1]
48real, parameter :: h43 = 369.4389437509002*(i_cp0*i_ts**3) ! Tp to Tc fit zS**4 * T**3 coef. [degC-2]
49real, parameter :: h44 = -33.83664947895248*(i_cp0*i_ts**4) ! Tp to Tc fit zS**4 * T**4 coef. [degC-3]
50real, parameter :: h45 = -9.987880382780322*(i_cp0*i_ts**5) ! Tp to Tc fit zS**4 * T**5 coef. [degC-4]
51real, parameter :: h50 = 246.9598888781377*i_cp0 ! Tp to Tc fit zS**5 coef. [degC]
52real, parameter :: h60 = 123.59576582457964*i_cp0 ! Tp to Tc fit zS**6 coef. [degC]
53real, parameter :: h70 = -48.5891069025409*i_cp0 ! Tp to Tc fit zS**7 coef. [degC]
54
55!>@}
56
57contains
58
59!> Convert input potential temperature [degC] and absolute salinity [g kg-1] to returned
60!! conservative temperature [degC] using the polynomial expressions from TEOS-10.
61elemental real function potemp_to_constemp(T, Sa) result(Tc)
62 real, intent(in) :: t !< Potential temperature [degC]
63 real, intent(in) :: sa !< Absolute salinity [g kg-1]
64
65 ! Local variables
66 real :: x2 ! Absolute salinity normalized by a plausible salinity range [nondim]
67 real :: x ! Square root of normalized absolute salinity [nondim]
68
69 x2 = max(i_s0 * sa, 0.0)
70 x = sqrt(x2)
71
72 tc = h00 + (t*(h01 + t*(h02 + t*(h03 + t*(h04 + t*(h05 + t*(h06 + t* h07)))))) &
73 + x2*(h20 + (t*(h21 + t*(h22 + t*(h23 + t*(h24 + t*(h25 + t*h26))))) &
74 + x*(h30 + (t*(h31 + t*(h32 + t*(h33 + t* h34))) &
75 + x*(h40 + (t*(h41 + t*(h42 + t*(h43 + t*(h44 + t*h45)))) &
76 + x*(h50 + x*(h60 + x* h70)) )) )) )) )
77
78end function potemp_to_constemp
79
80
81!> Return the partial derivative of conservative temperature with potential temperature [nondim]
82!! based on the polynomial expressions from TEOS-10.
83elemental real function dtc_dtp(T, Sa)
84 real, intent(in) :: t !< Potential temperature [degC]
85 real, intent(in) :: sa !< Absolute salinity [g kg-1]
86
87 ! Local variables
88 real :: x2 ! Absolute salinity normalized by a plausible salinity range [nondim]
89 real :: x ! Square root of normalized absolute salinity [nondim]
90
91 x2 = max(i_s0 * sa, 0.0)
92 x = sqrt(x2)
93
94 dtc_dtp = ( h01 + t*(2.*h02 + t*(3.*h03 + t*(4.*h04 + t*(5.*h05 + t*(6.*h06 + t*(7.*h07)))))) ) &
95 + x2*( (h21 + t*(2.*h22 + t*(3.*h23 + t*(4.*h24 + t*(5.*h25 + t*(6.*h26)))))) &
96 + x*( (h31 + t*(2.*h32 + t*(3.*h33 + t*(4.*h34)))) &
97 + x*(h41 + t*(2.*h42 + t*(3.*h43 + t*(4.*h44 + t*(5.*h45))))) ) )
98
99end function dtc_dtp
100
101
102
103!> Convert input potential temperature [degC] and absolute salinity [g kg-1] to returned
104!! conservative temperature [degC] by inverting the polynomial expressions from TEOS-10.
105elemental real function constemp_to_potemp(Tc, Sa) result(Tp)
106 real, intent(in) :: tc !< Conservative temperature [degC]
107 real, intent(in) :: sa !< Absolute salinity [g kg-1]
108
109 real :: tp_num ! The numerator of a simple expression for potential temperature [degC]
110 real :: i_tp_den ! The inverse of the denominator of a simple expression for potential temperature [nondim]
111 real :: tc_diff ! The difference between an estimate of conservative temperature and its target [degC]
112 real :: tp_old ! A previous estimate of the potential tempearture [degC]
113 real :: dtp_dtc ! The partial derivative of potential temperature with conservative temperature [nondim]
114 ! The following are coefficients in the nominator (TPNxx) or denominator (TPDxx) of a simple rational
115 ! expression that approximately converts conservative temperature to potential temperature.
116 real, parameter :: tpn00 = -1.446013646344788e-2 ! Simple fit numerator constant [degC]
117 real, parameter :: tpn10 = -3.305308995852924e-3*sprac_sref ! Simple fit numerator Sa coef. [degC ppt-1]
118 real, parameter :: tpn20 = 1.062415929128982e-4*sprac_sref**2 ! Simple fit numerator Sa**2 coef. [degC ppt-2]
119 real, parameter :: tpn01 = 9.477566673794488e-1 ! Simple fit numerator Tc coef. [nondim]
120 real, parameter :: tpn11 = 2.166591947736613e-3*sprac_sref ! Simple fit numerator Sa * Tc coef. [ppt-1]
121 real, parameter :: tpn02 = 3.828842955039902e-3 ! Simple fit numerator Tc**2 coef. [degC-1]
122 real, parameter :: tpd10 = 6.506097115635800e-4*sprac_sref ! Simple fit denominator Sa coef. [ppt-1]
123 real, parameter :: tpd01 = 3.830289486850898e-3 ! Simple fit denominator Tc coef. [degC-1]
124 real, parameter :: tpd02 = 1.247811760368034e-6 ! Simple fit denominator Tc**2 coef. [degC-2]
125
126 ! Estimate the potential temperature and its derivative from an approximate rational function fit.
127 tp_num = tpn00 + (sa*(tpn10 + tpn20*sa) + tc*(tpn01 + (tpn11*sa + tpn02*tc)))
128 i_tp_den = 1.0 / (1.0 + (tpd10*sa + tc*(tpd01 + tpd02*tc)))
129 tp = tp_num*i_tp_den
130 dtp_dtc = ((tpn01 + (tpn11*sa + 2.*tpn02*tc)) - (tpd01 + 2.*tpd02*tc)*tp)*i_tp_den
131
132 ! Start the 1.5 iterations through the modified Newton-Raphson iterative method, which is also known
133 ! as the Newton-McDougall method. In this case 1.5 iterations converge to 64-bit machine precision
134 ! for oceanographically relevant temperatures and salinities.
135
136 tc_diff = potemp_to_constemp(tp, sa) - tc
137 tp_old = tp
138 tp = tp_old - tc_diff*dtp_dtc
139
140 dtp_dtc = 1.0 / dtc_dtp(0.5*(tp + tp_old), sa)
141
142 tp = tp_old - tc_diff*dtp_dtc
143 tc_diff = potemp_to_constemp(tp, sa) - tc
144 tp_old = tp
145
146 tp = tp_old - tc_diff*dtp_dtc
147
148end function constemp_to_potemp
149
150!> \namespace MOM_temperature_conv
151!!
152!! \section MOM_temperature_conv Temperature conversions
153!!
154!! This module has functions that convert potential temperature to conservative temperature
155!! and the reverse, as described in the TEOS-10 manual. This code was originally derived
156!! from their corresponding routines in the gsw code package, but has had some refactoring so that the
157!! answers are more likely to reproduce across compilers and levels of optimization. A complete
158!! discussion of the thermodynamics of seawater and the definition of conservative temperature
159!! can be found in IOC et al. (2010).
160!!
161!! \subsection section_temperature_conv_references References
162!!
163!! IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of seawater - 2010:
164!! Calculation and use of thermodynamic properties. Intergovernmental Oceanographic Commission,
165!! Manuals and Guides No. 56, UNESCO (English), 196 pp.
166!! (Available from www.teos-10.org/pubs/TEOS-10_Manual.pdf)
167