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

Last change on this file since 3538 was 3447, checked in by eckhard, 6 years ago

inifor: Renamed source files for compatibilty with PALM build system

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