source: palm/trunk/UTIL/inifor/src/util.f90 @ 3208

Last change on this file since 3208 was 3183, checked in by suehring, 6 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 10.4 KB
Line 
1!> @file src/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-2018 Leibniz Universitaet Hannover
18! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: util.f90 3183 2018-07-27 14:25:55Z sward $
28! Improved real-to-string conversion
29!
30!
31! 3182 2018-07-27 13:36:03Z suehring
32! Initial revision
33!
34!
35!
36! Authors:
37! --------
38! @author Eckhard Kadasch
39!
40! Description:
41! ------------
42!> The util module provides miscellaneous utility routines for INIFOR.
43!------------------------------------------------------------------------------!
44 MODULE util
45
46    USE, INTRINSIC :: ISO_C_BINDING,                                           &
47        ONLY :  C_CHAR, C_INT, C_PTR, C_SIZE_T
48    USE defs,                                                                  &
49        ONLY :  dp, PI, DATE, SNAME
50
51    IMPLICIT NONE
52
53    TYPE, BIND(c) :: tm_struct
54       INTEGER(C_INT) :: tm_sec     !< seconds after the minute [0, 61]
55       INTEGER(C_INT) :: tm_min     !< minutes after the hour [0, 59]
56       INTEGER(C_INT) :: tm_hour    !< hours since midnight [0, 23]
57       INTEGER(C_INT) :: tm_mday    !< day of the month [1, 31]
58       INTEGER(C_INT) :: tm_mon     !< month since January [0, 11]
59       INTEGER(C_INT) :: tm_year    !< years since 1900
60       INTEGER(C_INT) :: tm_wday    !< days since Sunday [0, 6]
61       INTEGER(C_INT) :: tm_yday    !< days since January 1st [0, 356]
62       INTEGER(C_INT) :: tm_isdst   !< Daylight Saving Time flag
63    END TYPE
64
65    INTERFACE
66
67       FUNCTION strptime(string, format, timeinfo) BIND(c, NAME='strptime')
68          IMPORT :: C_CHAR, C_SIZE_T, tm_struct
69
70          IMPLICIT NONE
71
72          CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN) ::  string, format
73          TYPE(tm_struct), INTENT(OUT)                     ::  timeinfo
74
75          INTEGER(C_SIZE_T)                                ::  strptime
76       END FUNCTION
77
78
79       FUNCTION strftime(string, string_len, format, timeinfo) BIND(c, NAME='strftime')
80          IMPORT :: C_CHAR, C_SIZE_T, tm_struct
81
82          IMPLICIT NONE
83
84          CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(OUT) ::  string
85          CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN)  ::  format
86          INTEGER(C_SIZE_T), INTENT(IN)                     ::  string_len
87          TYPE(tm_struct), INTENT(IN)                       ::  timeinfo
88
89          INTEGER(C_SIZE_T)                                 ::  strftime
90       END FUNCTION
91
92
93       FUNCTION mktime(timeinfo) BIND(c, NAME='mktime')
94          IMPORT :: C_PTR, tm_struct
95
96          IMPLICIT NONE
97
98          TYPE(tm_struct), INTENT(IN) ::  timeinfo
99
100          TYPE(C_PTR)                 ::  mktime
101       END FUNCTION
102
103    END INTERFACE
104
105 CONTAINS
106
107    CHARACTER(LEN=DATE) FUNCTION add_hours_to(date_string, hours)
108       CHARACTER(LEN=DATE), INTENT(IN)          ::  date_string
109       INTEGER, INTENT(IN)                      ::  hours
110
111       CHARACTER(KIND=C_CHAR, LEN=*), PARAMETER ::  format_string = "%Y%m%d%H"
112       CHARACTER(KIND=C_CHAR, LEN=DATE)         ::  c_date_string
113       TYPE(C_PTR)                              ::  c_pointer
114       TYPE(tm_struct)                          ::  time_info
115       INTEGER                                  ::  err
116
117       c_date_string = date_string
118
119       ! Convert C string to C tm struct
120       CALL init_tm(time_info)
121       err = strptime(c_date_string, format_string, time_info)
122   
123       ! Manipulate and normalize C tm struct
124       time_info % tm_hour = time_info % tm_hour + hours
125       c_pointer = mktime(time_info)
126
127       ! Convert back to C string
128       err = strftime(c_date_string, INT(DATE, KIND=C_SIZE_T),                 &
129                      format_string, time_info)
130
131       add_hours_to = c_date_string
132    END FUNCTION
133
134
135    SUBROUTINE print_tm(timeinfo)
136       TYPE(tm_struct), INTENT(IN) :: timeinfo
137
138       PRINT *, "sec: ", timeinfo % tm_sec,  &  !< seconds after the minute [0, 61]
139                "min: ", timeinfo % tm_min,  &  !< minutes after the hour [0, 59]
140                "hr:  ", timeinfo % tm_hour, &  !< hours since midnight [0, 23]
141                "day: ", timeinfo % tm_mday, &  !< day of the month [1, 31]
142                "mon: ", timeinfo % tm_mon,  &  !< month since January [0, 11]
143                "yr:  ", timeinfo % tm_year, &  !< years since 1900
144                "wday:", timeinfo % tm_wday, &  !< days since Sunday [0, 6]
145                "yday:", timeinfo % tm_yday, &  !< days since January 1st [0, 356]
146                "dst: ", timeinfo % tm_isdst    !< Daylight Saving time flag
147    END SUBROUTINE print_tm
148
149   
150    SUBROUTINE init_tm(timeinfo)
151       TYPE(tm_struct), INTENT(INOUT) :: timeinfo
152
153       timeinfo % tm_sec   = 0
154       timeinfo % tm_min   = 0
155       timeinfo % tm_hour  = 0
156       timeinfo % tm_mday  = 0
157       timeinfo % tm_mon   = 0
158       timeinfo % tm_year  = 0
159       timeinfo % tm_wday  = 0
160       timeinfo % tm_yday  = 0
161
162       ! We use UTC times, so marking Daylight Saving Time (DST) 'not available'
163       ! (< 0). If this is set to 0, mktime will convert the timeinfo to DST and
164       ! add one hour.
165       timeinfo % tm_isdst = -1
166    END SUBROUTINE init_tm
167
168
169    SUBROUTINE fake_output_3d(a)
170
171       REAL(dp), INTENT(INOUT)       ::  a(:,:,:)
172       REAL(dp)                      ::  lxi, lyi
173       INTEGER ::  i,j,k
174
175       lyi = 2.0_dp * PI / (SIZE(a, 2) - 1.0_dp)
176       lxi = 2.0_dp * PI / (SIZE(a, 1) - 1.0_dp)
177
178       DO k = 1, SIZE(a, 3)
179       DO j = 1, SIZE(a, 2)
180       DO i = 1, SIZE(a, 1)
181           a(i,j,k) = SIN(lxi * i) * COS(lyi * j) + k
182       END DO
183       END DO
184       END DO
185
186    END SUBROUTINE fake_output_3d
187
188
189    SUBROUTINE fake_output_2d(a, offset)
190
191       REAL(dp), INTENT(INOUT) ::  a(:,:)
192       INTEGER, INTENT(IN)     ::  offset
193       REAL(dp)                ::  lxi, lyi
194       INTEGER                 ::  i,j
195
196       lyi = 2.0_dp*PI / (SIZE(a, 2) - 1.0_dp)
197       lxi = 2.0_dp*PI / (SIZE(a, 1) - 1.0_dp)
198
199       a(:,:) = 1.0_dp
200       DO j = 1, SIZE(a, 2)
201       DO i = 1, SIZE(a, 1)
202          a(i,j) = SIN(lxi * i) * COS(lyi * j) + offset
203       END DO
204       END DO
205
206    END SUBROUTINE fake_output_2d
207
208
209    SUBROUTINE linspace(start, stop, array)
210
211       REAL(dp), INTENT(IN)    ::  start, stop
212       REAL(dp), INTENT(INOUT) ::  array(0:)
213       INTEGER                 ::  i, n
214
215       n = UBOUND(array, 1)
216
217       IF (n .EQ. 0)  THEN
218
219          array(0) = start
220
221       ELSE
222
223          DO i = 0, n
224             array(i) = start + REAL(i, dp) / n * (stop - start)
225          END DO
226
227       END IF
228       
229    END SUBROUTINE linspace
230
231
232!------------------------------------------------------------------------------!
233! Description:
234! ------------
235!> Reverse the order of the third (vertical) array dimension from top-down
236!> (COSMO) to bottom-up (PALM)
237!------------------------------------------------------------------------------!
238    SUBROUTINE reverse(input_arr)
239
240       REAL(dp), INTENT(INOUT) ::  input_arr(:,:,:)
241
242       input_arr = input_arr(:,:,size(input_arr, 3):1:-1)
243
244    END SUBROUTINE reverse
245
246
247    SUBROUTINE deaverage(avg_1, t1, avg_2, t2, avg_3, t3)
248
249       REAL(dp), DIMENSION(:,:,:), INTENT(IN)  ::  avg_1, avg_2
250       REAL(dp), INTENT(IN)                    ::  t1, t2, t3
251       REAL(dp), DIMENSION(:,:,:), INTENT(OUT) ::  avg_3
252
253       REAL(dp)                                ::  ti
254 
255       ti = 1.0_dp / t3
256
257       avg_3(:,:,:) = ti * ( t2 * avg_2(:,:,:) - t1 * avg_1(:,:,:) )
258
259    END SUBROUTINE deaverage
260
261
262    SUBROUTINE get_basic_state(z, beta, p_sl, t_sl, rd, g, p0)
263
264       REAL(dp), INTENT(IN)  ::  z(1:)  !< height [m]
265       REAL(dp), INTENT(IN)  ::  beta   !< logarithmic lapse rate, dT / d ln(p) [K]
266       REAL(dp), INTENT(IN)  ::  p_sl   !< reference pressure [Pa]
267       REAL(dp), INTENT(IN)  ::  t_sl   !< reference tempereature [K]
268       REAL(dp), INTENT(IN)  ::  rd     !< ideal gas constant of dry air [J/kg/K]
269       REAL(dp), INTENT(IN)  ::  g      !< acceleration of Earth's gravity [m/s^2]
270       REAL(dp), INTENT(OUT) ::  p0(1:) !< COSMO-DE basic state pressure [Pa]
271       REAL(dp) ::  root_frac, factor   !< precomputed factors
272
273       factor = - t_sl / beta
274       root_frac = (2.0_dp * beta * g) / (rd * t_sl*t_sl)
275
276       p0(:) = p_sl * EXP(                                                     &
277                  factor * ( 1.0_dp - SQRT( 1.0_dp - root_frac * z(:) ) )  &
278               )
279
280    END SUBROUTINE get_basic_state
281
282
283    ! Convert a real number to a string in scientific notation
284    ! showing four significant digits.
285    CHARACTER(LEN=SNAME) FUNCTION real_to_str(val, format)
286
287        REAL(dp), INTENT(IN)                   ::  val
288        CHARACTER(LEN=*), OPTIONAL, INTENT(IN) ::  format
289
290        IF (PRESENT(format))  THEN
291           WRITE(real_to_str, format) val
292        ELSE
293           WRITE(real_to_str, '(E11.4)') val
294        END IF
295        real_to_str = ADJUSTL(real_to_str)
296
297    END FUNCTION real_to_str
298
299
300    CHARACTER(LEN=12) FUNCTION real_to_str_f(val)
301
302        REAL(dp), INTENT(IN) ::  val
303
304        WRITE(real_to_str_f, '(F12.4)') val
305        real_to_str_f = ADJUSTL(real_to_str_f)
306
307    END FUNCTION real_to_str_f
308
309
310    CHARACTER(LEN=10) FUNCTION str(val)
311
312        INTEGER, INTENT(IN) ::  val
313
314        WRITE(str, '(i10)') val
315        str = ADJUSTL(str)
316
317    END FUNCTION str
318
319
320    CHARACTER(LEN=30) FUNCTION strs(vals)
321
322        INTEGER, INTENT(IN) ::  vals(:)
323        INTEGER ::  i
324
325        strs = ''
326        DO i = 1, SIZE(vals)
327           strs = TRIM(strs) // TRIM(str(vals(i)))
328        END DO
329
330    END FUNCTION strs
331
332
333    SUBROUTINE normalize_path(path)
334       
335        CHARACTER(LEN=*), INTENT(INOUT) ::  path
336        INTEGER ::  n
337
338        n = LEN_TRIM(path)
339
340        IF (path(n:n) .NE. '/')  THEN
341           path = TRIM(path) // '/'
342        END IF
343
344    END SUBROUTINE
345
346 END MODULE
347
Note: See TracBrowser for help on using the repository browser.