source: palm/trunk/UTIL/inifor/src/inifor_util.f90

Last change on this file was 4843, checked in by raasch, 21 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

  • Property svn:keywords set to Id
File size: 16.9 KB
Line 
1!> @file src/inifor_util.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2021 Leibniz Universitaet Hannover
18! Copyright 2017-2021 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor_util.f90 4843 2021-01-15 15:22:11Z banzhafs $
28! Support integer-string conversion for 8 and 4 bit types
29!
30!
31! 4568 2020-06-19 11:56:30Z eckhard
32! Added function for checking floating point equality
33!
34!
35! 4523 2020-05-07 15:58:16Z eckhard
36! respect integer working precision (iwp) specified in inifor_defs.f90
37!
38!
39! 4499 2020-04-16 15:51:56Z eckhard
40! Bugfix: avoid using already overwritten elements in 'reverse' subroutine
41!
42!
43! 4481 2020-03-31 18:55:54Z maronga
44! Bugfix: use explicit loop in 'reverse' subroutine instead of implicit loop
45!
46! 3866 2019-04-05 14:25:01Z eckhard
47! Use PALM's working precision
48! Improved coding style
49!
50!
51! 3785 2019-03-06 10:41:14Z eckhard
52! Prefixed all INIFOR modules with inifor_
53!
54!
55! 3557 2018-11-22 16:01:22Z eckhard
56! Updated documentation
57!
58!
59! 3447 2018-10-29 15:52:54Z eckhard
60! Renamed source files for compatibilty with PALM build system
61!
62!
63! 3395 2018-10-22 17:32:49Z eckhard
64! New routines for computing potential temperature and moist air density
65! Increased number of digits in real-to-str conversion
66!
67! 3183 2018-07-27 14:25:55Z suehring
68! Improved real-to-string conversion
69!
70!
71! 3182 2018-07-27 13:36:03Z suehring
72! Initial revision
73!
74!
75!
76! Authors:
77! --------
78!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
79!
80! Description:
81! ------------
82!> The util module provides miscellaneous utility routines for INIFOR.
83!------------------------------------------------------------------------------!
84 MODULE inifor_util
85
86    USE inifor_defs,                                                           &
87        ONLY :  PI, DATE, SNAME, iwp, wp
88    USE inifor_types,                                                          &
89        ONLY :  grid_definition
90    USE, INTRINSIC :: ISO_C_BINDING,                                           &
91        ONLY :  C_CHAR, C_INT, C_PTR, C_SIZE_T
92
93    IMPLICIT NONE
94
95!------------------------------------------------------------------------------!
96! Description:
97! ------------
98!> Fortran implementation of C's struct tm for representing points in time
99!------------------------------------------------------------------------------!
100    TYPE, BIND(c) :: tm_struct
101       INTEGER(C_INT) :: tm_sec     !< seconds after the minute [0, 61]
102       INTEGER(C_INT) :: tm_min     !< minutes after the hour [0, 59]
103       INTEGER(C_INT) :: tm_hour    !< hours since midnight [0, 23]
104       INTEGER(C_INT) :: tm_mday    !< day of the month [1, 31]
105       INTEGER(C_INT) :: tm_mon     !< month since January [0, 11]
106       INTEGER(C_INT) :: tm_year    !< years since 1900
107       INTEGER(C_INT) :: tm_wday    !< days since Sunday [0, 6]
108       INTEGER(C_INT) :: tm_yday    !< days since January 1st [0, 356]
109       INTEGER(C_INT) :: tm_isdst   !< Daylight Saving Time flag
110    END TYPE
111
112 INTERFACE
113
114!------------------------------------------------------------------------------!
115! Description:
116! ------------
117!> Interface to C's strptime function, which converts a string to a tm
118!> structure.
119!------------------------------------------------------------------------------!
120    FUNCTION strptime(string, format, timeinfo) BIND(c, NAME='strptime')
121       IMPORT :: C_CHAR, C_SIZE_T, tm_struct
122
123       IMPLICIT NONE
124
125       CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN) ::  string, format
126       TYPE(tm_struct), INTENT(OUT)                     ::  timeinfo
127
128       INTEGER(C_SIZE_T)                                ::  strptime
129    END FUNCTION
130
131
132!------------------------------------------------------------------------------!
133! Description:
134! ------------
135!> Interface to C's strftime function, which converts the given 'timeinfo'
136!> structure to a string in the given 'format'.
137!------------------------------------------------------------------------------!
138    FUNCTION strftime(string, string_len, format, timeinfo) BIND(c, NAME='strftime')
139       IMPORT :: C_CHAR, C_SIZE_T, tm_struct
140
141       IMPLICIT NONE
142
143       CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(OUT) ::  string
144       CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN)  ::  format
145       INTEGER(C_SIZE_T), INTENT(IN)                     ::  string_len
146       TYPE(tm_struct), INTENT(IN)                       ::  timeinfo
147
148       INTEGER(C_SIZE_T)                                 ::  strftime
149    END FUNCTION
150
151
152!------------------------------------------------------------------------------!
153! Description:
154! ------------
155!> Interface to C's mktime function, which converts the given tm structure to
156!> Unix time (number of seconds since 0 UTC, 1st January 1970). INIFOR uses the
157!> side effect that a call to mktime noramlizes the given 'timeinfo' structure,
158!> e.g. increments the date if hours overfow 24.
159!------------------------------------------------------------------------------!
160    FUNCTION mktime(timeinfo) BIND(c, NAME='mktime')
161       IMPORT :: C_PTR, tm_struct
162
163       IMPLICIT NONE
164
165       TYPE(tm_struct), INTENT(IN) ::  timeinfo
166
167       TYPE(C_PTR)                 ::  mktime
168    END FUNCTION
169
170 END INTERFACE
171
172 
173 INTERFACE str
174    MODULE PROCEDURE str_iwp
175    MODULE PROCEDURE str_4
176 END INTERFACE
177
178
179 CONTAINS
180
181!------------------------------------------------------------------------------!
182! Description:
183! ------------
184!> Takes a string of the form YYYYMMDDHH, adds the given number of hours to its
185!> tm structure representation, and returns the corresponding string in the same
186!> format
187!------------------------------------------------------------------------------!
188 CHARACTER(LEN=DATE) FUNCTION add_hours_to(date_string, hours)
189    CHARACTER(LEN=DATE), INTENT(IN)          ::  date_string
190    INTEGER(iwp), INTENT(IN)                 ::  hours
191
192    CHARACTER(KIND=C_CHAR, LEN=*), PARAMETER ::  format_string = "%Y%m%d%H"
193    CHARACTER(KIND=C_CHAR, LEN=DATE)         ::  c_date_string
194    TYPE(C_PTR)                              ::  c_pointer
195    TYPE(tm_struct)                          ::  time_info
196    INTEGER(iwp)                             ::  err
197
198    c_date_string = date_string
199
200    ! Convert C string to C tm struct
201    CALL init_tm(time_info)
202    err = strptime(c_date_string, format_string, time_info)
203 
204    ! Manipulate and normalize C tm struct
205    time_info%tm_hour = time_info%tm_hour + hours
206    c_pointer = mktime(time_info)
207
208    ! Convert back to C string
209    err = strftime(c_date_string, INT(DATE, KIND=C_SIZE_T),                 &
210                   format_string, time_info)
211
212    add_hours_to = c_date_string
213 END FUNCTION
214
215
216!------------------------------------------------------------------------------!
217! Description:
218! ------------
219!> Print all members of the given tm structure
220!------------------------------------------------------------------------------!
221 SUBROUTINE print_tm(timeinfo)
222    TYPE(tm_struct), INTENT(IN) :: timeinfo
223
224    PRINT *, "sec: ", timeinfo%tm_sec,  &  !< seconds after the minute [0, 61]
225             "min: ", timeinfo%tm_min,  &  !< minutes after the hour [0, 59]
226             "hr:  ", timeinfo%tm_hour, &  !< hours since midnight [0, 23]
227             "day: ", timeinfo%tm_mday, &  !< day of the month [1, 31]
228             "mon: ", timeinfo%tm_mon,  &  !< month since January [0, 11]
229             "yr:  ", timeinfo%tm_year, &  !< years since 1900
230             "wday:", timeinfo%tm_wday, &  !< days since Sunday [0, 6]
231             "yday:", timeinfo%tm_yday, &  !< days since January 1st [0, 356]
232             "dst: ", timeinfo%tm_isdst    !< Daylight Saving time flag
233 END SUBROUTINE print_tm
234
235   
236!------------------------------------------------------------------------------!
237! Description:
238! ------------
239!> Initialize the given tm structure with zero values
240!------------------------------------------------------------------------------!
241 SUBROUTINE init_tm(timeinfo)
242    TYPE(tm_struct), INTENT(INOUT) :: timeinfo
243
244    timeinfo%tm_sec   = 0
245    timeinfo%tm_min   = 0
246    timeinfo%tm_hour  = 0
247    timeinfo%tm_mday  = 0
248    timeinfo%tm_mon   = 0
249    timeinfo%tm_year  = 0
250    timeinfo%tm_wday  = 0
251    timeinfo%tm_yday  = 0
252
253    ! We use UTC times, so marking Daylight Saving Time (DST) 'not available'
254    ! (< 0). If this is set to 0, mktime will convert the timeinfo to DST and
255    ! add one hour.
256    timeinfo%tm_isdst = -1
257 END SUBROUTINE init_tm
258
259
260!------------------------------------------------------------------------------!
261! Description:
262! ------------
263!> Fill the given array with values equally spaced between and including start
264!> and stop
265!------------------------------------------------------------------------------!
266 SUBROUTINE linspace(start, stop, array)
267
268    REAL(wp), INTENT(IN)    ::  start, stop
269    REAL(wp), INTENT(INOUT) ::  array(0:)
270    INTEGER(iwp)            ::  i, n
271
272    n = UBOUND(array, 1)
273
274    IF (n .EQ. 0)  THEN
275
276       array(0) = start
277
278    ELSE
279
280       DO i = 0, n
281          array(i) = start + REAL(i, wp) / n * (stop - start)
282       ENDDO
283
284    ENDIF
285   
286 END SUBROUTINE linspace
287
288
289!------------------------------------------------------------------------------!
290! Description:
291! ------------
292!> Reverse the order of the third (vertical) array dimension from top-down
293!> (COSMO) to bottom-up (PALM)
294!------------------------------------------------------------------------------!
295 SUBROUTINE reverse(input_arr)
296
297    INTEGER(iwp) ::  idx, opposite_idx, half_idx
298    INTEGER(iwp) ::  size_1st_dimension
299    INTEGER(iwp) ::  size_2nd_dimension
300    INTEGER(iwp) ::  size_3rd_dimension
301    INTEGER(iwp) ::  lbound_3rd_dimension
302    INTEGER(iwp) ::  ubound_3rd_dimension
303
304    REAL(wp), INTENT(INOUT) ::  input_arr(:,:,:)
305    REAL(wp), ALLOCATABLE  :: buffer_arr(:,:)
306
307    lbound_3rd_dimension = LBOUND(input_arr, 3)
308    ubound_3rd_dimension = UBOUND(input_arr, 3)
309    size_1st_dimension = SIZE(input_arr, 1)
310    size_2nd_dimension = SIZE(input_arr, 2)
311    size_3rd_dimension = SIZE(input_arr, 3)
312    half_idx = lbound_3rd_dimension + size_3rd_dimension / 2 - 1
313
314    ALLOCATE( buffer_arr(size_1st_dimension, size_2nd_dimension) )
315
316    DO  idx = lbound_3rd_dimension, half_idx
317       opposite_idx = ubound_3rd_dimension - ( idx - lbound_3rd_dimension )
318       buffer_arr(:,:) = input_arr(:,:,idx)
319       input_arr(:,:,idx) = input_arr(:,:,opposite_idx)
320       input_arr(:,:,opposite_idx) = buffer_arr(:,:)
321    ENDDO
322   
323    DEALLOCATE( buffer_arr )
324
325 END SUBROUTINE reverse
326
327
328!------------------------------------------------------------------------------!
329! Description:
330! ------------
331!>
332!------------------------------------------------------------------------------!
333 SUBROUTINE deaverage(avg_1, t1, avg_2, t2, avg_3, t3)
334
335    REAL(wp), DIMENSION(:,:,:), INTENT(IN)  ::  avg_1, avg_2
336    REAL(wp), INTENT(IN)                    ::  t1, t2, t3
337    REAL(wp), DIMENSION(:,:,:), INTENT(OUT) ::  avg_3
338
339    REAL(wp)                                ::  ti
340 
341    ti = 1.0_wp / t3
342
343    avg_3(:,:,:) = ti * ( t2 * avg_2(:,:,:) - t1 * avg_1(:,:,:) )
344
345 END SUBROUTINE deaverage
346
347
348!------------------------------------------------------------------------------!
349! Description:
350! ------------
351!> Compute the COSMO-DE/-D2 basic state pressure profile
352!------------------------------------------------------------------------------!
353 SUBROUTINE get_basic_state(z, beta, p_sl, t_sl, rd, g, p0)
354
355    REAL(wp), INTENT(IN)  ::  z(1:)  !< height [m]
356    REAL(wp), INTENT(IN)  ::  beta   !< logarithmic lapse rate, dT / d ln(p) [K]
357    REAL(wp), INTENT(IN)  ::  p_sl   !< reference pressure [Pa]
358    REAL(wp), INTENT(IN)  ::  t_sl   !< reference tempereature [K]
359    REAL(wp), INTENT(IN)  ::  rd     !< ideal gas constant of dry air [J/kg/K]
360    REAL(wp), INTENT(IN)  ::  g      !< acceleration of Earth's gravity [m/s^2]
361    REAL(wp), INTENT(OUT) ::  p0(1:) !< COSMO basic state pressure [Pa]
362    REAL(wp) ::  root_frac, factor   !< precomputed factors
363
364    factor = - t_sl / beta
365    root_frac = (2.0_wp * beta * g) / (rd * t_sl*t_sl)
366
367    p0(:) = p_sl * EXP(                                                     &
368               factor * ( 1.0_wp - SQRT( 1.0_wp - root_frac * z(:) ) )      &
369    )
370
371 END SUBROUTINE get_basic_state
372
373
374!------------------------------------------------------------------------------!
375! Description:
376! ------------
377!> Converts the absolute temperature to the potential temperature in place using
378!> the identity a^b = e^(b ln(a)).
379!>
380!>     theta = T * (p_ref/p)^(R/c_p) = T * e^( R/c_p * ln(p_ref/p) )
381!------------------------------------------------------------------------------!
382 SUBROUTINE potential_temperature(t, p, p_ref, r, cp)
383    REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  t
384    REAL(wp), DIMENSION(:,:,:), INTENT(IN)    ::  p
385    REAL(wp), INTENT(IN)                      ::  p_ref, r, cp
386    REAL(wp)                                  ::  rcp
387
388    rcp = r/cp
389    t(:,:,:) =  t(:,:,:) * EXP( rcp * LOG(p_ref / p(:,:,:)) )
390
391 END SUBROUTINE potential_temperature
392
393
394!------------------------------------------------------------------------------!
395! Description:
396! ------------
397!> Compute the density in place of the given temperature (t_rho).
398!------------------------------------------------------------------------------!
399 SUBROUTINE moist_density(t_rho, p, qv, rd, rv)
400    REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  t_rho
401    REAL(wp), DIMENSION(:,:,:), INTENT(IN)    ::  p, qv
402    REAL(wp), INTENT(IN)                      ::  rd, rv
403
404    t_rho(:,:,:) = p(:,:,:) / (                                             &
405       (rv * qv(:,:,:) + rd * (1.0_wp - qv(:,:,:))) * t_rho(:,:,:)          &
406    )
407
408 END SUBROUTINE moist_density
409
410!------------------------------------------------------------------------------!
411! Description:
412! ------------
413! Convert a real number to a string in scientific notation showing four
414! significant digits.
415!------------------------------------------------------------------------------!
416 CHARACTER(LEN=SNAME) FUNCTION real_to_str(val, format)
417
418    REAL(wp), INTENT(IN)                   ::  val
419    CHARACTER(LEN=*), OPTIONAL, INTENT(IN) ::  format
420
421    IF (PRESENT( format ) )  THEN
422       WRITE( real_to_str, format ) val
423    ELSE
424       WRITE( real_to_str, '(E11.4)' ) val
425    ENDIF
426    real_to_str = ADJUSTL( real_to_str )
427
428 END FUNCTION real_to_str
429
430
431!------------------------------------------------------------------------------!
432! Description:
433! ------------
434!> Converts the given real value to a string
435!------------------------------------------------------------------------------!
436 CHARACTER(LEN=16) FUNCTION real_to_str_f(val)
437
438     REAL(wp), INTENT(IN) ::  val
439
440     WRITE(real_to_str_f, '(F16.8)') val
441     real_to_str_f = ADJUSTL(real_to_str_f)
442
443 END FUNCTION real_to_str_f
444
445
446!------------------------------------------------------------------------------!
447! Description:
448! ------------
449!> Converts the given integer value to a string
450!------------------------------------------------------------------------------!
451 CHARACTER(LEN=10) FUNCTION str_iwp(val)
452
453     INTEGER(iwp), INTENT(IN) ::  val
454
455     WRITE(str_iwp, '(i10)') val
456     str_iwp= ADJUSTL(str_iwp)
457
458 END FUNCTION str_iwp
459
460
461 CHARACTER(LEN=10) FUNCTION str_4(val)
462
463     INTEGER(4), INTENT(IN) ::  val
464
465     WRITE(str_4, '(i10)') val
466     str_4 = ADJUSTL(str_4)
467
468 END FUNCTION str_4
469!------------------------------------------------------------------------------!
470! Description:
471! ------------
472!> If the given path is not conlcuded by a slash, add one.
473!------------------------------------------------------------------------------!
474 SUBROUTINE normalize_path(path)
475     
476     CHARACTER(LEN=*), INTENT(INOUT) ::  path
477     INTEGER(iwp) ::  n
478
479     n = LEN_TRIM(path)
480
481     IF (path(n:n) .NE. '/')  THEN
482        path = TRIM(path) // '/'
483     ENDIF
484
485 END SUBROUTINE
486
487
488!------------------------------------------------------------------------------!
489! Description:
490! ------------
491!> Check for exact or near equality of two floating point numbers. Inspired by
492!> https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
493!------------------------------------------------------------------------------!
494 LOGICAL FUNCTION nearly_equal(a, b, max_rel_diff)
495
496    REAL(wp), INTENT(IN) ::  a, b, max_rel_diff
497    REAL(wp)             ::  diff, mag
498 
499    diff = ABS( a - b ) 
500    mag = MAX( ABS(a), ABS(b) )
501    nearly_equal = ( diff .LE. mag * max_rel_diff )
502
503 END FUNCTION nearly_equal
504
505 END MODULE inifor_util
Note: See TracBrowser for help on using the repository browser.