SICOPOLIS V5-dev  Revision 1105
enth_temp_omega_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : e n t h _ t e m p _ o m e g a _ m
4 !
5 !> @file
6 !!
7 !! Conversion from temperature (temp) and water content (omega) to enthalpy
8 !! (enth) and vice versa.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2013-2017 Ralf Greve, Heinz Blatter
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Conversion from temperature (temp) and water content (omega) to enthalpy
35 !! (enth) and vice versa.
36 !<------------------------------------------------------------------------------
38 
39 use sico_types_m
40 
41 implicit none
42 save
43 
44 !> c_int_table: Temperature integral of the specific heat of ice.
45 !> Index is temperature in deg C.
46  real(dp), dimension(-256:255), private :: c_int_table
47 
48 !> c_int_inv_table: Inverse of the temperature integral of the specific heat
49 !> of ice. Index is enthalpy in J/kg (zero for 0 deg C).
50  real(dp), dimension(-524288:524287), private :: c_int_inv_table
51 
52 !> n_temp_min: Lower index limit of properly defined values in c_int_table
53 !> (n_temp_min >= -256).
54  integer(i4b), private :: n_temp_min
55 
56 !> n_temp_max: Upper index limit of properly defined values in c_int_table
57 !> (n_temp_max <= 255).
58  integer(i4b), private :: n_temp_max
59 
60 !> n_enth_min: Lower index limit of properly defined values in c_int_inv_table
61 !> (n_enth_min >= -524288).
62  integer(i4b), private :: n_enth_min
63 
64 !> n_enth_max: Upper index limit of properly defined values in c_int_inv_table
65 !> (n_enth_max <= 524287).
66  integer(i4b), private :: n_enth_max
67 
68 !> L: Latent heat of ice
69  real(dp), private :: l
70 
71 !> L_inv: Inverse of the latent heat of ice
72  real(dp), private :: l_inv
73 
76 private :: c_int_val, c_int_inv_val
77 
78 contains
79 
80 !-------------------------------------------------------------------------------
81 !> Computation of the temperature integral of the specific heat of ice as a
82 !! table (c_int_table). Further, definition of the latent heat of ice.
83 !<------------------------------------------------------------------------------
84 subroutine calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
85 
86 implicit none
87 
88 integer(i4b), intent(in) :: n_tmp_min, n_tmp_max
89 real(dp), dimension(n_tmp_min:n_tmp_max), intent(in) :: c_table
90 real(dp), optional, intent(in) :: L_val
91 
92 integer(i4b) :: n
93 real(dp) :: c_int_zero
94 
95 !-------- Initialisation --------
96 
97 c_int_table = 0.0_dp
98 
99 n_temp_min = n_tmp_min
100 n_temp_max = n_tmp_max
101 
102 if ((n_temp_min <= -256).or.(n_temp_max >= 255)) &
103  stop ' >>> calc_c_int_table: Temperature indices out of allowed range!'
104 
105 !-------- Numerical integration with the trapezoidal rule (spacing
106 ! of data in c_table and c_int_table assumed to be 1 deg C) --------
107 
108 do n=n_temp_min+1, n_temp_max
109  c_int_table(n) = c_int_table(n-1) + 0.5_dp*(c_table(n-1)+c_table(n))
110  ! that's the real stuff
111 end do
112 
113 do n=n_temp_max+1, 255
114  c_int_table(n) = c_int_table(n_temp_max) ! dummy values
115 end do
116 
117 !-------- Shift of the zero level to 0 deg C --------
118 
119 c_int_zero = c_int_table(0)
120 
121 do n=-256, 255
122  c_int_table(n) = c_int_table(n) - c_int_zero
123 end do
124 
125 !-------- Further stuff: Latent heat of ice --------
126 
127 if ( present(l_val) ) then
128  l = l_val
129 else
130  l = 3.35e+05_dp ! in J/kg
131 end if
132 
133 l_inv = 1.0_dp/l
134 
135 end subroutine calc_c_int_table
136 
137 !-------------------------------------------------------------------------------
138 !> Computation of the inverse of the temperature integral of the specific heat
139 !! of ice as a table (c_int_inv_table).
140 !<------------------------------------------------------------------------------
141 subroutine calc_c_int_inv_table()
143 implicit none
144 
145 integer(i4b) :: n
146 integer(i4b) :: n_temp_1, n_temp_2
147 real(dp) :: enth_min, enth_max
148 real(dp) :: enth_val, enth_1, enth_2
149 
150 !-------- Initialisation --------
151 
152 c_int_inv_table = 0.0_dp
153 
154 enth_min = c_int_val(real(n_temp_min,dp))
155 enth_max = c_int_val(real(n_temp_max,dp))
156 
157 n_enth_min = ceiling(enth_min)
158 n_enth_max = floor(enth_max)
159 
160 if ((n_enth_min <= -524288).or.(n_enth_max >= 524287)) &
161  stop ' >>> calc_c_int_inv_table: Enthalpy indices out of allowed range!'
162 
163 !-------- Linear interpolation between adjacent enthalpy values --------
164 
165 n_temp_1 = n_temp_min
166 n_temp_2 = n_temp_min+1
167 
168 do n=n_enth_min, n_enth_max
169 
170  enth_val = real(n,dp)
171 
172  do
173 
174  if ((n_temp_1 > n_temp_max).or.(n_temp_2 > n_temp_max)) &
175  stop ' >>> calc_c_int_inv_table: Temperature indices out of allowed range!'
176 
177  enth_1 = c_int_val(real(n_temp_1,dp))
178  enth_2 = c_int_val(real(n_temp_2,dp))
179 
180  if ( (enth_1 <= enth_val).and.(enth_2 >= enth_val) ) exit
181 
182  n_temp_1 = n_temp_1+1
183  n_temp_2 = n_temp_2+1
184 
185  end do
186 
187  c_int_inv_table(n) = real(n_temp_1,dp) &
188  + (real(n_temp_2,dp)-real(n_temp_1,dp)) &
189  * (enth_val-enth_1)/(enth_2-enth_1)
190  ! Linear interpolation
191 end do
192 
193 do n=-524288, n_enth_min-1
194  c_int_inv_table(n) = c_int_inv_table(n_enth_min) ! dummy values
195 end do
196 
197 do n=n_enth_max+1, 524287
198  c_int_inv_table(n) = c_int_inv_table(n_enth_max) ! dummy values
199 end do
200 
201 end subroutine calc_c_int_inv_table
202 
203 !-------------------------------------------------------------------------------
204 !> Temperature integral of the specific heat of ice
205 !! (enthalpy as function of temperature).
206 !<------------------------------------------------------------------------------
207 function c_int_val(temp_val)
208 
209 implicit none
210 
211 real(dp) :: c_int_val
212 
213 real(dp), intent(in) :: temp_val
214 
215 integer(i4b) :: n_temp_1, n_temp_2
216 
217 n_temp_1 = floor(temp_val)
218 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
219 n_temp_2 = n_temp_1 + 1
220 
221 c_int_val = c_int_table(n_temp_1) &
222  + (c_int_table(n_temp_2)-c_int_table(n_temp_1)) &
223  * (temp_val-real(n_temp_1,dp)) ! Linear interpolation
224 
225 end function c_int_val
226 
227 !-------------------------------------------------------------------------------
228 !> Inverse function of c_int_val (temperature as function of enthalpy).
229 !<------------------------------------------------------------------------------
230 function c_int_inv_val(enth_val)
231 
232 implicit none
233 
234 real(dp) :: c_int_inv_val
235 
236 real(dp), intent(in) :: enth_val
237 
238 integer(i4b) :: n_enth_1, n_enth_2
239 
240 n_enth_1 = floor(enth_val)
241 n_enth_1 = max(min(n_enth_1, n_enth_max-1), n_enth_min)
242 n_enth_2 = n_enth_1 + 1
243 
244 c_int_inv_val = c_int_inv_table(n_enth_1) &
245  + (c_int_inv_table(n_enth_2)-c_int_inv_table(n_enth_1)) &
246  * (enth_val-real(n_enth_1,dp)) ! Linear interpolation
247 
248 end function c_int_inv_val
249 
250 !-------------------------------------------------------------------------------
251 !> Enthalpy as a function of temperature and water content.
252 !<------------------------------------------------------------------------------
253 function enth_fct_temp_omega(temp_val, omega_val)
255 implicit none
256 
257 real(dp) :: enth_fct_temp_omega
258 
259 real(dp), intent(in) :: temp_val, omega_val
260 
261 enth_fct_temp_omega = c_int_val(temp_val) + l*omega_val
262 
263 end function enth_fct_temp_omega
264 
265 !-------------------------------------------------------------------------------
266 !> Temperature as a function of enthalpy.
267 !<------------------------------------------------------------------------------
268 function temp_fct_enth(enth_val, temp_m_val)
270 implicit none
271 
272 real(dp) :: temp_fct_enth
273 
274 real(dp), intent(in) :: enth_val
275 real(dp), intent(in) :: temp_m_val
276 
277 real(dp) :: enth_i
278 
279 enth_i = c_int_val(temp_m_val) ! Enthalpy of pure ice at the melting point
280 
281 if (enth_val < enth_i) then ! cold ice
282  temp_fct_enth = c_int_inv_val(enth_val)
283 else ! temperate ice
284  temp_fct_enth = temp_m_val
285 end if
286 
287 end function temp_fct_enth
288 
289 !-------------------------------------------------------------------------------
290 !> Water content as a function of enthalpy.
291 !<------------------------------------------------------------------------------
292 function omega_fct_enth(enth_val, temp_m_val)
294 implicit none
295 
296 real(dp) :: omega_fct_enth
297 
298 real(dp), intent(in) :: enth_val
299 real(dp), intent(in) :: temp_m_val
300 
301 real(dp) :: enth_i
302 
303 enth_i = c_int_val(temp_m_val) ! Enthalpy of pure ice at the melting point
304 
305 omega_fct_enth = max((enth_val-enth_i)*l_inv, 0.0_dp)
306 
307 end function omega_fct_enth
308 
309 !-------------------------------------------------------------------------------
310 
311 end module enth_temp_omega_m
312 !
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
real(dp) function, public temp_fct_enth(enth_val, temp_m_val)
Temperature as a function of enthalpy.
real(dp) function, public omega_fct_enth(enth_val, temp_m_val)
Water content as a function of enthalpy.
Declarations of kind types for SICOPOLIS.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
integer, parameter dp
Double-precision reals.
subroutine, public calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
Computation of the temperature integral of the specific heat of ice as a table (c_int_table). Further, definition of the latent heat of ice.