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

Last change on this file since 4004 was 3866, checked in by eckhard, 6 years ago

inifor: Use PALM's working precision; improved error handling, coding style, and comments

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