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

Last change on this file since 2698 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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