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

Last change on this file since 4499 was 4499, checked in by eckhard, 3 years ago

bugfix for explicit loop in 'reverse' subroutine, updated test suite

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