SICOPOLIS V5-dev  Revision 1288
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-2018 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 use error_m
41 
42 implicit none
43 save
44 
45 #ifndef ALLOW_OPENAD /* Normal */
46 
47 !> c_int_table: Temperature integral of the specific heat of ice.
48 !> Index is temperature in deg C.
49  real(dp), dimension(-256:255), private :: c_int_table
50 
51 !> c_int_inv_table: Inverse of the temperature integral of the specific heat
52 !> of ice. Index is enthalpy in J/kg (zero for 0 deg C).
53  real(dp), dimension(-524288:524287), private :: c_int_inv_table
54 
55 !> n_temp_min: Lower index limit of properly defined values in c_int_table
56 !> (n_temp_min >= -256).
57  integer(i4b), private :: n_temp_min
58 
59 !> n_temp_max: Upper index limit of properly defined values in c_int_table
60 !> (n_temp_max <= 255).
61  integer(i4b), private :: n_temp_max
62 
63 !> n_enth_min: Lower index limit of properly defined values in c_int_inv_table
64 !> (n_enth_min >= -524288).
65  integer(i4b), private :: n_enth_min
66 
67 !> n_enth_max: Upper index limit of properly defined values in c_int_inv_table
68 !> (n_enth_max <= 524287).
69  integer(i4b), private :: n_enth_max
70 
71 !> L: Latent heat of ice
72  real(dp), private :: l
73 
74 !> L_inv: Inverse of the latent heat of ice
75  real(dp), private :: l_inv
76 
79 private :: c_int_val, c_int_inv_val
80 
81 #else /* OpenAD */
82 
83 real(dp), dimension(-256:255), public :: c_int_table
84 real(dp), dimension(-524288:524287), public :: c_int_inv_table
85 integer(i4b), public :: n_temp_min
86 integer(i4b), public :: n_temp_max
87 integer(i4b), public :: n_enth_min
88 integer(i4b), public :: n_enth_max
89 real(dp), public :: l_inv
90 real(dp), public :: l_eto
91 
94 public :: c_int_val, c_int_inv_val
95 private :: myfloor_local_eto
96 
97 #endif /* Normal vs. OpenAD */
98 
99 contains
100 
101 !-------------------------------------------------------------------------------
102 !> Computation of the temperature integral of the specific heat of ice as a
103 !! table (c_int_table). Further, definition of the latent heat of ice.
104 !<------------------------------------------------------------------------------
105 subroutine calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
107 implicit none
108 
109 integer(i4b), intent(in) :: n_tmp_min, n_tmp_max
110 real(dp), dimension(n_tmp_min:n_tmp_max), intent(in) :: c_table
111 real(dp), optional, intent(in) :: L_val
112 
113 integer(i4b) :: n
114 real(dp) :: c_int_zero
115 character(len=256) :: errormsgg
116 
117 !-------- Initialisation --------
118 
119 c_int_table = 0.0_dp
120 
121 n_temp_min = n_tmp_min
122 n_temp_max = n_tmp_max
123 
124 if ((n_temp_min <= -256).or.(n_temp_max >= 255)) then
125  errormsgg = ' >>> calc_c_int_table: ' &
126  //'Temperature indices out of allowed range!'
127  call error(errormsgg)
128 end if
129 
130 !-------- Numerical integration with the trapezoidal rule (spacing
131 ! of data in c_table and c_int_table assumed to be 1 deg C) --------
132 
133 do n=n_temp_min+1, n_temp_max
134  c_int_table(n) = c_int_table(n-1) + 0.5_dp*(c_table(n-1)+c_table(n))
135  ! that's the real stuff
136 end do
137 
138 do n=n_temp_max+1, 255
139  c_int_table(n) = c_int_table(n_temp_max) ! dummy values
140 end do
141 
142 !-------- Shift of the zero level to 0 deg C --------
143 
144 c_int_zero = c_int_table(0)
145 
146 do n=-256, 255
147  c_int_table(n) = c_int_table(n) - c_int_zero
148 end do
149 
150 !-------- Further stuff: Latent heat of ice --------
151 
152 if ( present(l_val) ) then
153 
154 #ifndef ALLOW_OPENAD /* Normal */
155  l = l_val
156 #else /* OpenAD */
157  l_eto = l_val
158 #endif /* Normal vs. OpenAD */
159 
160 else
161 
162 #ifndef ALLOW_OPENAD /* Normal */
163  l = 3.35e+05_dp ! in J/kg
164 #else /* OpenAD */
165  l_eto = 3.35e+05_dp ! in J/kg
166 #endif /* Normal vs. OpenAD */
167 
168 end if
169 
170 #ifndef ALLOW_OPENAD /* Normal */
171  l_inv = 1.0_dp/l
172 #else /* OpenAD */
173  l_inv = 1.0_dp/l_eto
174 #endif /* Normal vs. OpenAD */
175 
176 end subroutine calc_c_int_table
177 
178 !-------------------------------------------------------------------------------
179 !> Computation of the inverse of the temperature integral of the specific heat
180 !! of ice as a table (c_int_inv_table).
181 !<------------------------------------------------------------------------------
182 subroutine calc_c_int_inv_table()
184 #ifdef ALLOW_OPENAD /* Normal */
185  use ctrl_m, only: myceiling, myfloor
186 #endif /* OpenAD */
187 
188 implicit none
189 
190 integer(i4b) :: n
191 integer(i4b) :: n_temp_1, n_temp_2
192 real(dp) :: enth_min, enth_max
193 real(dp) :: enth_val, enth_1, enth_2
194 character(len=256) :: errormsgg
195 
196 !-------- Initialisation --------
197 
198 c_int_inv_table = 0.0_dp
199 
200 enth_min = c_int_val(real(n_temp_min,dp))
201 enth_max = c_int_val(real(n_temp_max,dp))
202 
203 #ifndef ALLOW_OPENAD /* Normal */
204 n_enth_min = ceiling(enth_min)
205 #else /* OpenAD */
206 call myceiling(enth_min,n_enth_min)
207 #endif /* Normal vs. OpenAD */
208 
209 #ifndef ALLOW_OPENAD /* Normal */
210 n_enth_max = floor(enth_max)
211 #else /* OpenAD */
212 call myfloor(enth_max, n_enth_max)
213 #endif /* Normal vs. OpenAD */
214 
215 if ((n_enth_min <= -524288).or.(n_enth_max >= 524287)) then
216  errormsgg = ' >>> calc_c_int_inv_table: ' &
217  //'Enthalpy indices out of allowed range!'
218  call error(errormsgg)
219 end if
220 
221 !-------- Linear interpolation between adjacent enthalpy values --------
222 
223 n_temp_1 = n_temp_min
224 n_temp_2 = n_temp_min+1
225 
226 do n=n_enth_min, n_enth_max
227 
228  enth_val = real(n,dp)
229 
230 #ifndef ALLOW_OPENAD /* Normal */
231 
232  do
233 
234  if ((n_temp_1 > n_temp_max).or.(n_temp_2 > n_temp_max)) then
235  errormsgg = ' >>> calc_c_int_inv_table: ' &
236  //'Temperature indices out of allowed range!'
237  call error(errormsgg)
238  end if
239 
240  enth_1 = c_int_val(real(n_temp_1,dp))
241  enth_2 = c_int_val(real(n_temp_2,dp))
242 
243  if ( (enth_1 <= enth_val).and.(enth_2 >= enth_val) ) exit
244 
245  n_temp_1 = n_temp_1+1
246  n_temp_2 = n_temp_2+1
247 
248  end do
249 
250 #else /* OpenAD */
251 
252  enth_1 = c_int_val(real(n_temp_1,dp))
253  enth_2 = c_int_val(real(n_temp_2,dp))
254 
255  do while ( .not.((enth_1 <= enth_val).and.(enth_2 >= enth_val)) )
256 
257  if ((n_temp_1 > n_temp_max).or.(n_temp_2 > n_temp_max)) then
258  errormsgg = ' >>> calc_c_int_inv_table: ' &
259  //'Temperature indices out of allowed range!'
260  call error(errormsgg)
261  end if
262 
263  enth_1 = c_int_val(real(n_temp_1,dp))
264  enth_2 = c_int_val(real(n_temp_2,dp))
265 
266  if ( .not.((enth_1 <= enth_val).and.(enth_2 >= enth_val)) ) then
267  n_temp_1 = n_temp_1+1
268  n_temp_2 = n_temp_2+1
269  end if
270 
271  end do
272 
273 #endif /* Normal vs. OpenAD */
274 
275  c_int_inv_table(n) = real(n_temp_1,dp) &
276  + (real(n_temp_2,dp)-real(n_temp_1,dp)) &
277  * (enth_val-enth_1)/(enth_2-enth_1)
278  ! Linear interpolation
279 end do
280 
281 do n=-524288, n_enth_min-1
282  c_int_inv_table(n) = c_int_inv_table(n_enth_min) ! dummy values
283 end do
284 
285 do n=n_enth_max+1, 524287
286  c_int_inv_table(n) = c_int_inv_table(n_enth_max) ! dummy values
287 end do
288 
289 end subroutine calc_c_int_inv_table
290 
291 !-------------------------------------------------------------------------------
292 !> Temperature integral of the specific heat of ice
293 !! (enthalpy as function of temperature).
294 !<------------------------------------------------------------------------------
295 function c_int_val(temp_val)
296 
297 implicit none
298 
299 real(dp) :: c_int_val
300 
301 real(dp), intent(in) :: temp_val
302 
303 integer(i4b) :: n_temp_1, n_temp_2
304 
305 ! character(len=256) :: errormsgg
306 
307 #ifndef ALLOW_OPENAD /* Normal */
308 n_temp_1 = floor(temp_val)
309 #else /* OpenAD */
310 n_temp_1 = myfloor_local_eto(temp_val)
311 #endif /* Normal vs. OpenAD */
312 
313 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
314 n_temp_2 = n_temp_1 + 1
315 
316 ! if ((n_temp_1 < n_temp_min-1).or.(n_temp_2 > n_temp_max+1)) then
317 ! errormsgg = ' >>> c_int_val: Temperature argument out of allowed range!'
318 ! call error(errormsgg)
319 ! end if
320 ! *** Commented out after some testing in order to save computing time. ***
321 
322 c_int_val = c_int_table(n_temp_1) &
323  + (c_int_table(n_temp_2)-c_int_table(n_temp_1)) &
324  * (temp_val-real(n_temp_1,dp)) ! Linear interpolation
325 
326 end function c_int_val
327 
328 !-------------------------------------------------------------------------------
329 !> Inverse function of c_int_val (temperature as function of enthalpy).
330 !<------------------------------------------------------------------------------
331 function c_int_inv_val(enth_val)
332 
333 implicit none
334 
335 real(dp) :: c_int_inv_val
336 
337 real(dp), intent(in) :: enth_val
338 
339 integer(i4b) :: n_enth_1, n_enth_2
340 
341 ! character(len=256) :: errormsgg
342 
343 #ifndef ALLOW_OPENAD /* Normal */
344 n_enth_1 = floor(enth_val)
345 #else /* OpenAD */
346 n_enth_1 = myfloor_local_eto(enth_val)
347 #endif /* Normal vs. OpenAD */
348 
349 n_enth_1 = max(min(n_enth_1, n_enth_max-1), n_enth_min)
350 n_enth_2 = n_enth_1 + 1
351 
352 ! if ((n_enth_1 < n_enth_min-1).or.(n_enth_2 > n_enth_max+1)) then
353 ! errormsgg = ' >>> c_int_inv_val: Enthalpy argument out of allowed range!'
354 ! call error(errormsgg)
355 ! end if
356 ! *** Commented out after some testing in order to save computing time. ***
357 
358 c_int_inv_val = c_int_inv_table(n_enth_1) &
359  + (c_int_inv_table(n_enth_2)-c_int_inv_table(n_enth_1)) &
360  * (enth_val-real(n_enth_1,dp)) ! Linear interpolation
361 
362 end function c_int_inv_val
363 
364 !-------------------------------------------------------------------------------
365 !> Enthalpy as a function of temperature and water content.
366 !<------------------------------------------------------------------------------
367 function enth_fct_temp_omega(temp_val, omega_val)
369 implicit none
370 
371 real(dp) :: enth_fct_temp_omega
372 
373 real(dp), intent(in) :: temp_val, omega_val
374 
375 #ifndef ALLOW_OPENAD /* Normal */
376 enth_fct_temp_omega = c_int_val(temp_val) + l*omega_val
377 #else /* OpenAD */
378 enth_fct_temp_omega = c_int_val(temp_val) + l_eto*omega_val
379 #endif /* Normal vs. OpenAD */
380 
381 end function enth_fct_temp_omega
382 
383 !-------------------------------------------------------------------------------
384 !> Temperature as a function of enthalpy.
385 !<------------------------------------------------------------------------------
386 function temp_fct_enth(enth_val, temp_m_val)
388 implicit none
389 
390 real(dp) :: temp_fct_enth
391 
392 real(dp), intent(in) :: enth_val
393 real(dp), intent(in) :: temp_m_val
394 
395 real(dp) :: enth_i
396 
397 enth_i = c_int_val(temp_m_val) ! Enthalpy of pure ice at the melting point
398 
399 if (enth_val < enth_i) then ! cold ice
400  temp_fct_enth = c_int_inv_val(enth_val)
401 else ! temperate ice
402  temp_fct_enth = temp_m_val
403 end if
404 
405 end function temp_fct_enth
406 
407 !-------------------------------------------------------------------------------
408 !> Water content as a function of enthalpy.
409 !<------------------------------------------------------------------------------
410 function omega_fct_enth(enth_val, temp_m_val)
412 implicit none
413 
414 real(dp) :: omega_fct_enth
415 
416 real(dp), intent(in) :: enth_val
417 real(dp), intent(in) :: temp_m_val
418 
419 real(dp) :: enth_i
420 
421 enth_i = c_int_val(temp_m_val) ! Enthalpy of pure ice at the melting point
422 
423 omega_fct_enth = max((enth_val-enth_i)*l_inv, 0.0_dp)
424 
425 end function omega_fct_enth
426 
427 !-------------------------------------------------------------------------------
428 
429 #ifdef ALLOW_OPENAD /* OpenAD */
430 
431  function myfloor_local_eto(num)
432 
433  use sico_types_m
434  use sico_variables_m
435 
436  implicit none
437 
438  integer(i4b) :: inum
439  real(dp), intent(in) :: num
440  integer(i4b) :: myfloor_local_ETO
441 
442  inum = int(num);
443 
444  if (abs(num-real(inum,dp)) <= &
445  (abs(num+real(inum,dp))*myepsilon_dp) ) then
446  myfloor_local_eto = inum;
447  else if (num>0) then
448  myfloor_local_eto = inum;
449  else if (num<0) then
450  myfloor_local_eto = inum - 1;
451  end if
452 
453  end function myfloor_local_eto
454 
455 #endif /* OpenAD */
456 
457 !-------------------------------------------------------------------------------
458 
459 end module enth_temp_omega_m
460 !
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_...
Definition: ctrl_m.F90:9
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
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.
Writing of error messages and stopping execution.
Definition: error_m.F90:35
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.
Declarations of global variables for SICOPOLIS.