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

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

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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