Ignore:
Timestamp:
Oct 22, 2018 5:32:49 PM (6 years ago)
Author:
eckhard
Message:

inifor: Added computation of geostrophic winds from COSMO input

File:
1 edited

Legend:

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

    r3262 r3395  
    2626! -----------------
    2727! $Id$
    28 ! Removed unnecessary check for output file
     28! Added command-line options for configuring the computation of geostrophic
     29!     winds (--averagin-mode, --averaging-angle)
     30! Added command-line option --input-prefix for setting input file prefixes all
     31!     at once
     32! Added --debug option for more verbose terminal output
     33! Added option-specific *_is_set LOGICALs to indicate invocation from the
     34!     command-line
     35! Improved error messages in case of empty file-name strings
     36! Improved routine naming
    2937!
    3038! 3183 2018-07-27 14:25:55Z suehring
     
    148156       INTEGER                            ::  arg_count, i
    149157
     158       cfg % p0_is_set = .FALSE.
     159       cfg % ug_is_set = .FALSE.
     160       cfg % vg_is_set = .FALSE.
     161       cfg % flow_prefix_is_set = .FALSE.
     162       cfg % input_prefix_is_set = .FALSE.
     163       cfg % radiation_prefix_is_set = .FALSE.
     164       cfg % soil_prefix_is_set = .FALSE.
     165       cfg % soilmoisture_prefix_is_set = .FALSE.
     166
    150167       arg_count = COMMAND_ARGUMENT_COUNT()
    151168       IF (arg_count .GT. 0)  THEN
     
    170187             SELECT CASE( TRIM(option) )
    171188
     189             CASE( '--averaging-mode' )
     190                CALL get_option_argument( i, arg )
     191                cfg % averaging_mode = TRIM(arg)
     192
    172193             CASE( '-date', '-d', '--date' )
    173194                CALL get_option_argument( i, arg )
    174195                cfg % start_date = TRIM(arg)
     196
     197             CASE( '--debug' )
     198                cfg % debug = .TRUE.
    175199
    176200             ! Elevation of the PALM-4U domain above sea level
     
    181205             ! surface pressure, at z0
    182206             CASE( '-p0', '-r', '--surface-pressure' )
     207                cfg % p0_is_set = .TRUE.
    183208                CALL get_option_argument( i, arg )
    184209                READ(arg, *) cfg % p0
     
    186211             ! geostrophic wind in x direction
    187212             CASE( '-ug', '-u', '--geostrophic-u' )
     213                cfg % ug_is_set = .TRUE.
    188214                CALL get_option_argument( i, arg )
    189215                READ(arg, *) cfg % ug
     
    191217             ! geostrophic wind in y direction
    192218             CASE( '-vg', '-v', '--geostrophic-v' )
     219                cfg % vg_is_set = .TRUE.
    193220                CALL get_option_argument( i, arg )
    194221                READ(arg, *) cfg % vg
     
    208235             CASE( '-hhl', '-l', '--hhl-file' )
    209236                CALL get_option_argument( i, arg )
    210                  cfg % hhl_file = TRIM(arg)
     237                cfg % hhl_file = TRIM(arg)
     238
     239             CASE( '--input-prefix')
     240                CALL get_option_argument( i, arg )
     241                cfg % input_prefix = TRIM(arg)
     242                cfg % input_prefix_is_set = .TRUE.
     243   
     244             CASE( '-a', '--averaging-angle' )
     245                CALL get_option_argument( i, arg )
     246                READ(arg, *) cfg % averaging_angle
    211247
    212248             CASE( '-static', '-t', '--static-driver' )
    213249                CALL get_option_argument( i, arg )
    214                  cfg % static_driver_file = TRIM(arg)
     250                cfg % static_driver_file = TRIM(arg)
    215251
    216252             CASE( '-soil', '-s', '--soil-file')
    217253                CALL get_option_argument( i, arg )
    218                  cfg % soiltyp_file = TRIM(arg)
     254                cfg % soiltyp_file = TRIM(arg)
    219255
    220256             CASE( '--flow-prefix')
    221257                CALL get_option_argument( i, arg )
    222                  cfg % flow_prefix = TRIM(arg)
    223 
     258                cfg % flow_prefix = TRIM(arg)
     259                cfg % flow_prefix_is_set = .TRUE.
     260   
    224261             CASE( '--radiation-prefix')
    225262                CALL get_option_argument( i, arg )
    226                  cfg % radiation_prefix = TRIM(arg)
    227 
     263                cfg % radiation_prefix = TRIM(arg)
     264                cfg % radiation_prefix_is_set = .TRUE.
     265   
    228266             CASE( '--soil-prefix')
    229267                CALL get_option_argument( i, arg )
    230                  cfg % soil_prefix = TRIM(arg)
    231 
     268                cfg % soil_prefix = TRIM(arg)
     269                cfg % soil_prefix_is_set = .TRUE.
     270   
    232271             CASE( '--soilmoisture-prefix')
    233272                CALL get_option_argument( i, arg )
    234                  cfg % soilmoisture_prefix = TRIM(arg)
     273                cfg % soilmoisture_prefix = TRIM(arg)
     274                cfg % soilmoisture_prefix_is_set = .TRUE.
    235275
    236276             CASE( '-o', '--output' )
     
    297337
    298338      all_files_present = .TRUE.
    299       all_files_present = all_files_present .AND. file_present(cfg % hhl_file)
    300       all_files_present = all_files_present .AND. file_present(cfg % namelist_file)
    301       all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file)
     339      all_files_present = all_files_present .AND. file_present(cfg % hhl_file, 'HHL')
     340      all_files_present = all_files_present .AND. file_present(cfg % namelist_file, 'NAMELIST')
     341      all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file, 'SOILTYP')
    302342
    303343      ! Only check optional static driver file name, if it has been given.
    304344      IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
    305          all_files_present = all_files_present .AND. file_present(cfg % static_driver_file)
     345         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver')
    306346      END IF
    307347
     
    335375      END SELECT
    336376
     377      SELECT CASE( TRIM(cfg % averaging_mode) )
     378      CASE( 'level', 'height')
     379      CASE DEFAULT
     380         message = "Averaging mode '" // TRIM(cfg % averaging_mode) //&
     381                   "' is not supported. " //&
     382                   "Please select either 'height' or 'level', " //&
     383                   "or omit the --averaging-mode option entirely, which corresponds "//&
     384                   "to the latter."
     385         CALL abort( 'validate_config', message )
     386      END SELECT
     387
     388      IF ( cfg % ug_is_set .NEQV. cfg % vg_is_set )  THEN
     389         message = "You specified only one component of the geostrophic " // &
     390                   "wind. Please specify either both or none."
     391         CALL abort( 'validate_config', message )
     392      END IF
    337393
    338394   END SUBROUTINE validate_config
    339395
    340396
    341    LOGICAL FUNCTION file_present(filename)
     397   LOGICAL FUNCTION file_present(filename, kind)
    342398      CHARACTER(LEN=PATH), INTENT(IN) ::  filename
    343 
    344       INQUIRE(FILE=filename, EXIST=file_present)
    345 
    346       IF (.NOT. file_present)  THEN
    347          message = "The given file '" // "' does not exist."
     399      CHARACTER(LEN=*), INTENT(IN)    ::  kind
     400
     401      IF (LEN(TRIM(filename))==0)  THEN
     402
     403         file_present = .FALSE.
     404         message = "No name was given for the " // TRIM(kind) // " file."
    348405         CALL report('file_present', message)
     406
     407      ELSE
     408
     409         INQUIRE(FILE=filename, EXIST=file_present)
     410
     411         IF (.NOT. file_present)  THEN
     412            message = "The given " // TRIM(kind) // " file '" //               &
     413                      TRIM(filename) // "' does not exist."
     414            CALL report('file_present', message)
     415         END IF
     416
    349417      END IF
    350418
     
    388456#endif
    389457
    390 !
    391 !------------------------------------------------------------------------------
    392 !- Section 1: Write global NetCDF attributes
     458!------------------------------------------------------------------------------
     459!- Section 1: Define NetCDF dimensions and coordinates
     460!------------------------------------------------------------------------------
     461       nt = SIZE(output_file % time)
     462       nx = palm_grid % nx
     463       ny = palm_grid % ny
     464       nz = palm_grid % nz
     465       z0 = palm_grid % z0
     466
     467
     468!
     469!------------------------------------------------------------------------------
     470!- Section 2: Write global NetCDF attributes
    393471!------------------------------------------------------------------------------
    394472       CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
     
    406484       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
    407485       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
     486       ! FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
     487       ! FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
     488       ! FIXME: Standard v1.9., origin_z)
     489       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
    408490       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
    409491       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
    410492
    411493!
    412 !------------------------------------------------------------------------------
    413 !- Section 2: Define NetCDF dimensions and coordinates
    414 !------------------------------------------------------------------------------
    415        nt = SIZE(output_file % time)
    416        nx = palm_grid % nx
    417        ny = palm_grid % ny
    418        nz = palm_grid % nz
    419        z0 = palm_grid % z0
    420 
    421494!
    422495!------------------------------------------------------------------------------
     
    618691       ELSE
    619692
    620           nbuffers = SIZE( group % in_var_list )
     693          ! Allocate one input buffer per input_variable. If more quantities
     694          ! have to be computed than input variables exist in this group,
     695          ! allocate more buffers. For instance, in the thermodynamics group,
     696          ! there are three input variabels (PP, T, Qv) and four quantities
     697          ! necessart (P, Theta, Rho, qv) for the corresponding output fields
     698          ! (p0, Theta, qv, ug, and vg)
     699          nbuffers = MAX( group % n_inputs, group % n_output_quantities )
    621700          ALLOCATE( buffer(nbuffers) )
    622701 CALL run_control('time', 'alloc')
    623702         
    624           DO ivar = 1, nbuffers
     703          ! Read in all input variables, leave extra buffers-if any-untouched.
     704          DO ivar = 1, group % n_inputs
    625705
    626706             input_var => group % in_var_list(ivar)
     
    628708             ! Check wheather P or PP is present in input file
    629709             IF (input_var % name == 'P')  THEN
    630                 input_var % name = TRIM( get_pressure_var(input_file) )
     710                input_var % name = TRIM( get_pressure_varname(input_file) )
    631711 CALL run_control('time', 'read')
    632712             END IF
     
    668748!> perturbation, 'PP', and returns the appropriate string.
    669749!------------------------------------------------------------------------------!
    670     CHARACTER(LEN=2) FUNCTION get_pressure_var(input_file) RESULT(var)
     750    CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
    671751       CHARACTER(LEN=*) ::  input_file
    672752       INTEGER          ::  ncid, varid
     
    692772       CALL check(nf90_close(ncid))
    693773
    694     END FUNCTION get_pressure_var
     774    END FUNCTION get_pressure_varname
    695775
    696776
     
    805885                                   start=start(1:ndim+1) ) )
    806886
    807        CASE ( 'constant scalar profile' )
     887       CASE ( 'constant scalar profile', 'geostrophic', 'internal profile' )
    808888
    809889          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
Note: See TracChangeset for help on using the changeset viewer.