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

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

inifor: Updated documentation

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