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

Last change on this file since 4475 was 4475, checked in by gronemeier, 18 months ago

bugfixes/changes to INIFOR:

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