Ignore:
Timestamp:
Apr 5, 2019 2:25:01 PM (3 years ago)
Author:
eckhard
Message:

inifor: Use PALM's working precision; improved error handling, coding style, and comments

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/UTIL/inifor/src/inifor_util.f90

    r3785 r3866  
    2626! -----------------
    2727! $Id$
     28! Use PALM's working precision
     29! Improved coding style
     30!
     31!
     32! 3785 2019-03-06 10:41:14Z eckhard
    2833! Prefixed all INIFOR modules with inifor_
    2934!
     
    5863!> The util module provides miscellaneous utility routines for INIFOR.
    5964!------------------------------------------------------------------------------!
    60 #if defined ( __netcdf )
    6165 MODULE inifor_util
    6266
    6367    USE inifor_defs,                                                           &
    64         ONLY :  dp, PI, DATE, SNAME
     68        ONLY :  PI, DATE, SNAME, wp
    6569    USE inifor_types,                                                          &
    6670        ONLY :  grid_definition
     
    8791    END TYPE
    8892
    89     INTERFACE
     93 INTERFACE
    9094
    9195!------------------------------------------------------------------------------!
     
    9599!> structure.
    96100!------------------------------------------------------------------------------!
    97        FUNCTION strptime(string, format, timeinfo) BIND(c, NAME='strptime')
    98           IMPORT :: C_CHAR, C_SIZE_T, tm_struct
    99 
    100           IMPLICIT NONE
    101 
    102           CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN) ::  string, format
    103           TYPE(tm_struct), INTENT(OUT)                     ::  timeinfo
    104 
    105           INTEGER(C_SIZE_T)                                ::  strptime
    106        END FUNCTION
     101    FUNCTION strptime(string, format, timeinfo) BIND(c, NAME='strptime')
     102       IMPORT :: C_CHAR, C_SIZE_T, tm_struct
     103
     104       IMPLICIT NONE
     105
     106       CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN) ::  string, format
     107       TYPE(tm_struct), INTENT(OUT)                     ::  timeinfo
     108
     109       INTEGER(C_SIZE_T)                                ::  strptime
     110    END FUNCTION
    107111
    108112
     
    113117!> structure to a string in the given 'format'.
    114118!------------------------------------------------------------------------------!
    115        FUNCTION strftime(string, string_len, format, timeinfo) BIND(c, NAME='strftime')
    116           IMPORT :: C_CHAR, C_SIZE_T, tm_struct
    117 
    118           IMPLICIT NONE
    119 
    120           CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(OUT) ::  string
    121           CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN)  ::  format
    122           INTEGER(C_SIZE_T), INTENT(IN)                     ::  string_len
    123           TYPE(tm_struct), INTENT(IN)                       ::  timeinfo
    124 
    125           INTEGER(C_SIZE_T)                                 ::  strftime
    126        END FUNCTION
     119    FUNCTION strftime(string, string_len, format, timeinfo) BIND(c, NAME='strftime')
     120       IMPORT :: C_CHAR, C_SIZE_T, tm_struct
     121
     122       IMPLICIT NONE
     123
     124       CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(OUT) ::  string
     125       CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN)  ::  format
     126       INTEGER(C_SIZE_T), INTENT(IN)                     ::  string_len
     127       TYPE(tm_struct), INTENT(IN)                       ::  timeinfo
     128
     129       INTEGER(C_SIZE_T)                                 ::  strftime
     130    END FUNCTION
    127131
    128132
     
    135139!> e.g. increments the date if hours overfow 24.
    136140!------------------------------------------------------------------------------!
    137        FUNCTION mktime(timeinfo) BIND(c, NAME='mktime')
    138           IMPORT :: C_PTR, tm_struct
    139 
    140           IMPLICIT NONE
    141 
    142           TYPE(tm_struct), INTENT(IN) ::  timeinfo
    143 
    144           TYPE(C_PTR)                 ::  mktime
    145        END FUNCTION
    146 
    147     END INTERFACE
     141    FUNCTION mktime(timeinfo) BIND(c, NAME='mktime')
     142       IMPORT :: C_PTR, tm_struct
     143
     144       IMPLICIT NONE
     145
     146       TYPE(tm_struct), INTENT(IN) ::  timeinfo
     147
     148       TYPE(C_PTR)                 ::  mktime
     149    END FUNCTION
     150
     151 END INTERFACE
    148152
    149153 CONTAINS
     
    156160!> format
    157161!------------------------------------------------------------------------------!
    158     CHARACTER(LEN=DATE) FUNCTION add_hours_to(date_string, hours)
    159        CHARACTER(LEN=DATE), INTENT(IN)          ::  date_string
    160        INTEGER, INTENT(IN)                      ::  hours
    161 
    162        CHARACTER(KIND=C_CHAR, LEN=*), PARAMETER ::  format_string = "%Y%m%d%H"
    163        CHARACTER(KIND=C_CHAR, LEN=DATE)         ::  c_date_string
    164        TYPE(C_PTR)                              ::  c_pointer
    165        TYPE(tm_struct)                          ::  time_info
    166        INTEGER                                  ::  err
    167 
    168        c_date_string = date_string
    169 
    170        ! Convert C string to C tm struct
    171        CALL init_tm(time_info)
    172        err = strptime(c_date_string, format_string, time_info)
     162 CHARACTER(LEN=DATE) FUNCTION add_hours_to(date_string, hours)
     163    CHARACTER(LEN=DATE), INTENT(IN)          ::  date_string
     164    INTEGER, INTENT(IN)                      ::  hours
     165
     166    CHARACTER(KIND=C_CHAR, LEN=*), PARAMETER ::  format_string = "%Y%m%d%H"
     167    CHARACTER(KIND=C_CHAR, LEN=DATE)         ::  c_date_string
     168    TYPE(C_PTR)                              ::  c_pointer
     169    TYPE(tm_struct)                          ::  time_info
     170    INTEGER                                  ::  err
     171
     172    c_date_string = date_string
     173
     174    ! Convert C string to C tm struct
     175    CALL init_tm(time_info)
     176    err = strptime(c_date_string, format_string, time_info)
     177 
     178    ! Manipulate and normalize C tm struct
     179    time_info%tm_hour = time_info%tm_hour + hours
     180    c_pointer = mktime(time_info)
     181
     182    ! Convert back to C string
     183    err = strftime(c_date_string, INT(DATE, KIND=C_SIZE_T),                 &
     184                   format_string, time_info)
     185
     186    add_hours_to = c_date_string
     187 END FUNCTION
     188
     189
     190!------------------------------------------------------------------------------!
     191! Description:
     192! ------------
     193!> Print all members of the given tm structure
     194!------------------------------------------------------------------------------!
     195 SUBROUTINE print_tm(timeinfo)
     196    TYPE(tm_struct), INTENT(IN) :: timeinfo
     197
     198    PRINT *, "sec: ", timeinfo%tm_sec,  &  !< seconds after the minute [0, 61]
     199             "min: ", timeinfo%tm_min,  &  !< minutes after the hour [0, 59]
     200             "hr:  ", timeinfo%tm_hour, &  !< hours since midnight [0, 23]
     201             "day: ", timeinfo%tm_mday, &  !< day of the month [1, 31]
     202             "mon: ", timeinfo%tm_mon,  &  !< month since January [0, 11]
     203             "yr:  ", timeinfo%tm_year, &  !< years since 1900
     204             "wday:", timeinfo%tm_wday, &  !< days since Sunday [0, 6]
     205             "yday:", timeinfo%tm_yday, &  !< days since January 1st [0, 356]
     206             "dst: ", timeinfo%tm_isdst    !< Daylight Saving time flag
     207 END SUBROUTINE print_tm
     208
    173209   
    174        ! Manipulate and normalize C tm struct
    175        time_info % tm_hour = time_info % tm_hour + hours
    176        c_pointer = mktime(time_info)
    177 
    178        ! Convert back to C string
    179        err = strftime(c_date_string, INT(DATE, KIND=C_SIZE_T),                 &
    180                       format_string, time_info)
    181 
    182        add_hours_to = c_date_string
    183     END FUNCTION
    184 
    185 
    186 !------------------------------------------------------------------------------!
    187 ! Description:
    188 ! ------------
    189 !> Print all members of the given tm structure
    190 !------------------------------------------------------------------------------!
    191     SUBROUTINE print_tm(timeinfo)
    192        TYPE(tm_struct), INTENT(IN) :: timeinfo
    193 
    194        PRINT *, "sec: ", timeinfo % tm_sec,  &  !< seconds after the minute [0, 61]
    195                 "min: ", timeinfo % tm_min,  &  !< minutes after the hour [0, 59]
    196                 "hr:  ", timeinfo % tm_hour, &  !< hours since midnight [0, 23]
    197                 "day: ", timeinfo % tm_mday, &  !< day of the month [1, 31]
    198                 "mon: ", timeinfo % tm_mon,  &  !< month since January [0, 11]
    199                 "yr:  ", timeinfo % tm_year, &  !< years since 1900
    200                 "wday:", timeinfo % tm_wday, &  !< days since Sunday [0, 6]
    201                 "yday:", timeinfo % tm_yday, &  !< days since January 1st [0, 356]
    202                 "dst: ", timeinfo % tm_isdst    !< Daylight Saving time flag
    203     END SUBROUTINE print_tm
    204 
    205    
    206210!------------------------------------------------------------------------------!
    207211! Description:
     
    209213!> Initialize the given tm structure with zero values
    210214!------------------------------------------------------------------------------!
    211     SUBROUTINE init_tm(timeinfo)
    212        TYPE(tm_struct), INTENT(INOUT) :: timeinfo
    213 
    214        timeinfo % tm_sec   = 0
    215        timeinfo % tm_min   = 0
    216        timeinfo % tm_hour  = 0
    217        timeinfo % tm_mday  = 0
    218        timeinfo % tm_mon   = 0
    219        timeinfo % tm_year  = 0
    220        timeinfo % tm_wday  = 0
    221        timeinfo % tm_yday  = 0
    222 
    223        ! We use UTC times, so marking Daylight Saving Time (DST) 'not available'
    224        ! (< 0). If this is set to 0, mktime will convert the timeinfo to DST and
    225        ! add one hour.
    226        timeinfo % tm_isdst = -1
    227     END SUBROUTINE init_tm
     215 SUBROUTINE init_tm(timeinfo)
     216    TYPE(tm_struct), INTENT(INOUT) :: timeinfo
     217
     218    timeinfo%tm_sec   = 0
     219    timeinfo%tm_min   = 0
     220    timeinfo%tm_hour  = 0
     221    timeinfo%tm_mday  = 0
     222    timeinfo%tm_mon   = 0
     223    timeinfo%tm_year  = 0
     224    timeinfo%tm_wday  = 0
     225    timeinfo%tm_yday  = 0
     226
     227    ! We use UTC times, so marking Daylight Saving Time (DST) 'not available'
     228    ! (< 0). If this is set to 0, mktime will convert the timeinfo to DST and
     229    ! add one hour.
     230    timeinfo%tm_isdst = -1
     231 END SUBROUTINE init_tm
    228232
    229233
     
    234238!> and stop
    235239!------------------------------------------------------------------------------!
    236     SUBROUTINE linspace(start, stop, array)
    237 
    238        REAL(dp), INTENT(IN)    ::  start, stop
    239        REAL(dp), INTENT(INOUT) ::  array(0:)
    240        INTEGER                 ::  i, n
    241 
    242        n = UBOUND(array, 1)
    243 
    244        IF (n .EQ. 0)  THEN
    245 
    246           array(0) = start
    247 
    248        ELSE
    249 
    250           DO i = 0, n
    251              array(i) = start + REAL(i, dp) / n * (stop - start)
    252           ENDDO
    253 
    254        ENDIF
    255        
    256     END SUBROUTINE linspace
     240 SUBROUTINE linspace(start, stop, array)
     241
     242    REAL(wp), INTENT(IN)    ::  start, stop
     243    REAL(wp), INTENT(INOUT) ::  array(0:)
     244    INTEGER                 ::  i, n
     245
     246    n = UBOUND(array, 1)
     247
     248    IF (n .EQ. 0)  THEN
     249
     250       array(0) = start
     251
     252    ELSE
     253
     254       DO i = 0, n
     255          array(i) = start + REAL(i, wp) / n * (stop - start)
     256       ENDDO
     257
     258    ENDIF
     259   
     260 END SUBROUTINE linspace
    257261
    258262
     
    263267!> (COSMO) to bottom-up (PALM)
    264268!------------------------------------------------------------------------------!
    265     SUBROUTINE reverse(input_arr)
    266 
    267        REAL(dp), INTENT(INOUT) ::  input_arr(:,:,:)
    268 
    269        input_arr = input_arr(:,:,size(input_arr, 3):1:-1)
    270 
    271     END SUBROUTINE reverse
     269 SUBROUTINE reverse(input_arr)
     270
     271    REAL(wp), INTENT(INOUT) ::  input_arr(:,:,:)
     272
     273    input_arr = input_arr(:,:,size(input_arr, 3):1:-1)
     274
     275 END SUBROUTINE reverse
    272276
    273277
     
    277281!>
    278282!------------------------------------------------------------------------------!
    279     SUBROUTINE deaverage(avg_1, t1, avg_2, t2, avg_3, t3)
    280 
    281        REAL(dp), DIMENSION(:,:,:), INTENT(IN)  ::  avg_1, avg_2
    282        REAL(dp), INTENT(IN)                    ::  t1, t2, t3
    283        REAL(dp), DIMENSION(:,:,:), INTENT(OUT) ::  avg_3
    284 
    285        REAL(dp)                                ::  ti
     283 SUBROUTINE deaverage(avg_1, t1, avg_2, t2, avg_3, t3)
     284
     285    REAL(wp), DIMENSION(:,:,:), INTENT(IN)  ::  avg_1, avg_2
     286    REAL(wp), INTENT(IN)                    ::  t1, t2, t3
     287    REAL(wp), DIMENSION(:,:,:), INTENT(OUT) ::  avg_3
     288
     289    REAL(wp)                                ::  ti
    286290 
    287        ti = 1.0_dp / t3
    288 
    289        avg_3(:,:,:) = ti * ( t2 * avg_2(:,:,:) - t1 * avg_1(:,:,:) )
    290 
    291     END SUBROUTINE deaverage
     291    ti = 1.0_wp / t3
     292
     293    avg_3(:,:,:) = ti * ( t2 * avg_2(:,:,:) - t1 * avg_1(:,:,:) )
     294
     295 END SUBROUTINE deaverage
    292296
    293297
     
    297301!> Compute the COSMO-DE/-D2 basic state pressure profile
    298302!------------------------------------------------------------------------------!
    299     SUBROUTINE get_basic_state(z, beta, p_sl, t_sl, rd, g, p0)
    300 
    301        REAL(dp), INTENT(IN)  ::  z(1:)  !< height [m]
    302        REAL(dp), INTENT(IN)  ::  beta   !< logarithmic lapse rate, dT / d ln(p) [K]
    303        REAL(dp), INTENT(IN)  ::  p_sl   !< reference pressure [Pa]
    304        REAL(dp), INTENT(IN)  ::  t_sl   !< reference tempereature [K]
    305        REAL(dp), INTENT(IN)  ::  rd     !< ideal gas constant of dry air [J/kg/K]
    306        REAL(dp), INTENT(IN)  ::  g      !< acceleration of Earth's gravity [m/s^2]
    307        REAL(dp), INTENT(OUT) ::  p0(1:) !< COSMO basic state pressure [Pa]
    308        REAL(dp) ::  root_frac, factor   !< precomputed factors
    309 
    310        factor = - t_sl / beta
    311        root_frac = (2.0_dp * beta * g) / (rd * t_sl*t_sl)
    312 
    313        p0(:) = p_sl * EXP(                                                     &
    314                   factor * ( 1.0_dp - SQRT( 1.0_dp - root_frac * z(:) ) )      &
    315        )
    316 
    317     END SUBROUTINE get_basic_state
     303 SUBROUTINE get_basic_state(z, beta, p_sl, t_sl, rd, g, p0)
     304
     305    REAL(wp), INTENT(IN)  ::  z(1:)  !< height [m]
     306    REAL(wp), INTENT(IN)  ::  beta   !< logarithmic lapse rate, dT / d ln(p) [K]
     307    REAL(wp), INTENT(IN)  ::  p_sl   !< reference pressure [Pa]
     308    REAL(wp), INTENT(IN)  ::  t_sl   !< reference tempereature [K]
     309    REAL(wp), INTENT(IN)  ::  rd     !< ideal gas constant of dry air [J/kg/K]
     310    REAL(wp), INTENT(IN)  ::  g      !< acceleration of Earth's gravity [m/s^2]
     311    REAL(wp), INTENT(OUT) ::  p0(1:) !< COSMO basic state pressure [Pa]
     312    REAL(wp) ::  root_frac, factor   !< precomputed factors
     313
     314    factor = - t_sl / beta
     315    root_frac = (2.0_wp * beta * g) / (rd * t_sl*t_sl)
     316
     317    p0(:) = p_sl * EXP(                                                     &
     318               factor * ( 1.0_wp - SQRT( 1.0_wp - root_frac * z(:) ) )      &
     319    )
     320
     321 END SUBROUTINE get_basic_state
    318322
    319323
     
    326330!>     theta = T * (p_ref/p)^(R/c_p) = T * e^( R/c_p * ln(p_ref/p) )
    327331!------------------------------------------------------------------------------!
    328     SUBROUTINE potential_temperature(t, p, p_ref, r, cp)
    329        REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) ::  t
    330        REAL(dp), DIMENSION(:,:,:), INTENT(IN)    ::  p
    331        REAL(dp), INTENT(IN)                      ::  p_ref, r, cp
    332        REAL(dp)                                  ::  rcp
    333 
    334        rcp = r/cp
    335        t(:,:,:) =  t(:,:,:) * EXP( rcp * LOG(p_ref / p(:,:,:)) )
    336 
    337     END SUBROUTINE potential_temperature
     332 SUBROUTINE potential_temperature(t, p, p_ref, r, cp)
     333    REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  t
     334    REAL(wp), DIMENSION(:,:,:), INTENT(IN)    ::  p
     335    REAL(wp), INTENT(IN)                      ::  p_ref, r, cp
     336    REAL(wp)                                  ::  rcp
     337
     338    rcp = r/cp
     339    t(:,:,:) =  t(:,:,:) * EXP( rcp * LOG(p_ref / p(:,:,:)) )
     340
     341 END SUBROUTINE potential_temperature
    338342
    339343
     
    343347!> Compute the density in place of the given temperature (t_rho).
    344348!------------------------------------------------------------------------------!
    345    SUBROUTINE moist_density(t_rho, p, qv, rd, rv)
    346        REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) ::  t_rho
    347        REAL(dp), DIMENSION(:,:,:), INTENT(IN)    ::  p, qv
    348        REAL(dp), INTENT(IN)                      ::  rd, rv
    349 
    350        t_rho(:,:,:) = p(:,:,:) / (                                             &
    351           (rv * qv(:,:,:) + rd * (1.0_dp - qv(:,:,:))) * t_rho(:,:,:)          &
    352        )
    353 
    354     END SUBROUTINE moist_density
    355 
    356 
    357     ! Convert a real number to a string in scientific notation
    358     ! showing four significant digits.
    359     CHARACTER(LEN=SNAME) FUNCTION real_to_str(val, format)
    360 
    361         REAL(dp), INTENT(IN)                   ::  val
    362         CHARACTER(LEN=*), OPTIONAL, INTENT(IN) ::  format
    363 
    364         IF (PRESENT(format))  THEN
    365            WRITE(real_to_str, format) val
    366         ELSE
    367            WRITE(real_to_str, '(E11.4)') val
    368         ENDIF
    369         real_to_str = ADJUSTL(real_to_str)
    370 
    371     END FUNCTION real_to_str
     349 SUBROUTINE moist_density(t_rho, p, qv, rd, rv)
     350    REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  t_rho
     351    REAL(wp), DIMENSION(:,:,:), INTENT(IN)    ::  p, qv
     352    REAL(wp), INTENT(IN)                      ::  rd, rv
     353
     354    t_rho(:,:,:) = p(:,:,:) / (                                             &
     355       (rv * qv(:,:,:) + rd * (1.0_wp - qv(:,:,:))) * t_rho(:,:,:)          &
     356    )
     357
     358 END SUBROUTINE moist_density
     359
     360!------------------------------------------------------------------------------!
     361! Description:
     362! ------------
     363! Convert a real number to a string in scientific notation showing four
     364! significant digits.
     365!------------------------------------------------------------------------------!
     366 CHARACTER(LEN=SNAME) FUNCTION real_to_str(val, format)
     367
     368    REAL(wp), INTENT(IN)                   ::  val
     369    CHARACTER(LEN=*), OPTIONAL, INTENT(IN) ::  format
     370
     371    IF (PRESENT( format ) )  THEN
     372       WRITE( real_to_str, format ) val
     373    ELSE
     374       WRITE( real_to_str, '(E11.4)' ) val
     375    ENDIF
     376    real_to_str = ADJUSTL( real_to_str )
     377
     378 END FUNCTION real_to_str
    372379
    373380
     
    377384!> Converts the given real value to a string
    378385!------------------------------------------------------------------------------!
    379     CHARACTER(LEN=16) FUNCTION real_to_str_f(val)
    380 
    381         REAL(dp), INTENT(IN) ::  val
    382 
    383         WRITE(real_to_str_f, '(F16.8)') val
    384         real_to_str_f = ADJUSTL(real_to_str_f)
    385 
    386     END FUNCTION real_to_str_f
     386 CHARACTER(LEN=16) FUNCTION real_to_str_f(val)
     387
     388     REAL(wp), INTENT(IN) ::  val
     389
     390     WRITE(real_to_str_f, '(F16.8)') val
     391     real_to_str_f = ADJUSTL(real_to_str_f)
     392
     393 END FUNCTION real_to_str_f
    387394
    388395
     
    392399!> Converts the given integer value to a string
    393400!------------------------------------------------------------------------------!
    394     CHARACTER(LEN=10) FUNCTION str(val)
    395 
    396         INTEGER, INTENT(IN) ::  val
    397 
    398         WRITE(str, '(i10)') val
    399         str = ADJUSTL(str)
    400 
    401     END FUNCTION str
     401 CHARACTER(LEN=10) FUNCTION str(val)
     402
     403     INTEGER, INTENT(IN) ::  val
     404
     405     WRITE(str, '(i10)') val
     406     str = ADJUSTL(str)
     407
     408 END FUNCTION str
    402409
    403410
     
    407414!> If the given path is not conlcuded by a slash, add one.
    408415!------------------------------------------------------------------------------!
    409     SUBROUTINE normalize_path(path)
    410         
    411         CHARACTER(LEN=*), INTENT(INOUT) ::  path
    412         INTEGER ::  n
    413 
    414         n = LEN_TRIM(path)
    415 
    416         IF (path(n:n) .NE. '/')  THEN
    417            path = TRIM(path) // '/'
    418         ENDIF
    419 
    420     END SUBROUTINE
     416 SUBROUTINE normalize_path(path)
     417     
     418     CHARACTER(LEN=*), INTENT(INOUT) ::  path
     419     INTEGER ::  n
     420
     421     n = LEN_TRIM(path)
     422
     423     IF (path(n:n) .NE. '/')  THEN
     424        path = TRIM(path) // '/'
     425     ENDIF
     426
     427 END SUBROUTINE
    421428
    422429 END MODULE inifor_util
    423 #endif
    424 
Note: See TracChangeset for help on using the changeset viewer.