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

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

inifor: Added computation of geostrophic winds from COSMO input

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