Ignore:
Timestamp:
Jul 27, 2018 1:36:03 PM (6 years ago)
Author:
suehring
Message:

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

Location:
palm/trunk/UTIL/inifor/src
Files:
1 added
8 edited

Legend:

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

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
     23! Added version and copyright output
    2324!
    2425!
     
    4243
    4344    USE defs,                                                                  &
    44         ONLY:  LNAME, dp
     45        ONLY:  LNAME, dp, VERSION, COPYRIGHT
    4546
    4647    USE util,                                                                  &
     
    9596
    9697    END SUBROUTINE abort
     98
     99
     100    SUBROUTINE print_version()
     101       PRINT *, "INIFOR " // VERSION
     102       PRINT *, COPYRIGHT
     103    END SUBROUTINE print_version
    97104
    98105
  • palm/trunk/UTIL/inifor/src/defs.f90

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
     23! Updated defaults for soil extrapolation steps and nudging time-scale
     24! Improved handling of the start date string
     25! Added gas constant for water vapor
     26! Bumped INIFOR version number
    2327!
    2428!
     
    6266 REAL(dp), PARAMETER ::  T_SL = 288.15_dp                     !< Reference temperature for computation of COSMO-DE's basic state pressure [K]
    6367 REAL(dp), PARAMETER ::  BETA = 42.0_dp                       !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic state pressure [K]
    64  REAL(dp), PARAMETER ::  RD   = 287.05_dp                     !< specific gar constant of dry air, used in computation of COSMO-DE's basic state [J/kg/K]
     68 REAL(dp), PARAMETER ::  RD   = 287.05_dp                     !< specific gas constant of dry air, used in computation of COSMO-DE's basic state [J/kg/K]
     69 REAL(dp), PARAMETER ::  RV   = 461.51_dp                     !< specific gas constant of water vapor [J/kg/K]
    6570 REAL(dp), PARAMETER ::  G    = 9.80665_dp                    !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s]
    6671 REAL(dp), PARAMETER ::  RHO_L = 1e3_dp                       !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3], in [kg/m^3]
     
    7277
    7378 ! INIFOR parameters
    74  INTEGER, PARAMETER          ::  FILL_ITERATIONS = 10         !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-]
    75  REAL(dp), PARAMETER         ::  FORCING_FREQ = 3600.0_dp     !< Reference pressure for potential temperature [Pa]
    76  CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.1.4'            !< path to script for generating input file names
     79 INTEGER, PARAMETER          ::  FILL_ITERATIONS = 5          !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-]
     80 INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
     81 REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
     82 CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.3.0'            !< INIFOR version number
     83 CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2018 Leibniz Universitaet Hannover' // &
     84     NEW_LINE(' ') // ' Copyright 2017-2018 Deutscher Wetterdienst Offenbach' !< Copyright notice
    7785
    7886 END MODULE defs
  • palm/trunk/UTIL/inifor/src/grid.f90

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
    23 !
     23! Introduced new PALM grid stretching
     24! Updated variable names and metadata for PIDS v1.9 compatibility
     25! Better compatibility with older Intel compilers:
     26! - avoiding implicit array allocation with new get_netcdf_variable()
     27!   subroutine instead of function
     28! Improved command line interface:
     29! - Produce forcing variables when '-mode profile' is being used
     30! - Renamend initial-condition mode variable 'mode' to 'ic_mode'
     31! - Improved handling of the start date string
     32! Removed unnecessary variables and routines
     33!
    2434!
    2535! Former revisions:
     
    4757    USE defs,                                                                  &
    4858        ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp,     &
    49                SNAME, LNAME, PATH, FORCING_FREQ, WATER_ID, FILL_ITERATIONS,    &
     59               SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS,    &
    5060               BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L
    5161    USE io,                                                                    &
    5262        ONLY:  get_netcdf_variable, get_netcdf_attribute,                      &
    53                parse_command_line_arguments
     63               parse_command_line_arguments, validate_config
    5464    USE netcdf,                                                                &
    5565        ONLY:  NF90_MAX_NAME, NF90_MAX_VAR_DIMS
     
    6777    SAVE
    6878   
    69     REAL(dp) ::  phi_equat = 0.0_dp    !< latitude of rotated equator of COSMO-DE grid [rad]
    70     REAL(dp) ::  phi_n     = 0.0_dp    !< latitude of rotated pole of COSMO-DE grid [rad]
    71     REAL(dp) ::  lambda_n  = 0.0_dp    !< longitude of rotaded pole of COSMO-DE grid [rad]
    72     REAL(dp) ::  phi_c     = 0.0_dp    !< rotated-grid latitude of the center of the PALM domain [rad]
    73     REAL(dp) ::  lambda_c  = 0.0_dp    !< rotated-grid longitude of the centre of the PALM domain [rad]
    74     REAL(dp) ::  phi_cn    = 0.0_dp    !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
    75     REAL(dp) ::  lambda_cn = 0.0_dp    !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
    76     REAL(dp) ::  gam       = 0.0_dp    !< angle for working around phirot2phi/rlarot2rla bug
    77     REAL(dp) ::  dx        = 0.0_dp    !< PALM-4U grid spacing in x direction [m]
    78     REAL(dp) ::  dy        = 0.0_dp    !< PALM-4U grid spacing in y direction [m]
    79     REAL(dp) ::  dz        = 0.0_dp    !< PALM-4U grid spacing in z direction [m]
    80     REAL(dp) ::  dxi       = 0.0_dp    !< inverse PALM-4U grid spacing in x direction [m^-1]
    81     REAL(dp) ::  dyi       = 0.0_dp    !< inverse PALM-4U grid spacing in y direction [m^-1]
    82     REAL(dp) ::  dzi       = 0.0_dp    !< inverse PALM-4U grid spacing in z direction [m^-1]
    83     REAL(dp) ::  lx        = 0.0_dp    !< PALM-4U domain size in x direction [m]
    84     REAL(dp) ::  ly        = 0.0_dp    !< PALM-4U domain size in y direction [m]
    85     REAL(dp) ::  lz        = 0.0_dp    !< PALM-4U domain size in z direction [m]
    86     REAL(dp) ::  ug        = 0.0_dp    !< geostrophic wind in x direction [m/s]
    87     REAL(dp) ::  vg        = 0.0_dp    !< geostrophic wind in y direction [m/s]
    88     REAL(dp) ::  p0        = 0.0_dp    !< PALM-4U surface pressure, at z0 [Pa]
    89     REAL(dp) ::  x0        = 0.0_dp    !< x coordinate of PALM-4U Earth tangent [m]
    90     REAL(dp) ::  y0        = 0.0_dp    !< y coordinate of PALM-4U Earth tangent [m]
    91     REAL(dp) ::  z0        = 0.0_dp    !< Elevation of the PALM-4U domain above sea level [m]
    92     REAL(dp) ::  lonmin    = 0.0_dp    !< Minimunm longitude of COSMO-DE's rotated-pole grid
    93     REAL(dp) ::  lonmax    = 0.0_dp    !< Maximum longitude of COSMO-DE's rotated-pole grid
    94     REAL(dp) ::  latmin    = 0.0_dp    !< Minimunm latitude of COSMO-DE's rotated-pole grid
    95     REAL(dp) ::  latmax    = 0.0_dp    !< Maximum latitude of COSMO-DE's rotated-pole grid
    96     REAL(dp) ::  latitude  = 0.0_dp    !< geograpohical latitude of the PALM-4U origin, from inipar namelist [deg]
    97     REAL(dp) ::  longitude = 0.0_dp    !< geograpohical longitude of the PALM-4U origin, from inipar namelist [deg]
    98     REAL(dp) ::  origin_lat= 0.0_dp    !< geograpohical latitude of the PALM-4U origin, from static driver netCDF file [deg]
    99     REAL(dp) ::  origin_lon= 0.0_dp    !< geograpohical longitude of the PALM-4U origin, from static driver netCDF file [deg]
    100     REAL(dp) ::  end_time  = 0.0_dp    !< PALM-4U simulation time [s]
     79    REAL(dp) ::  phi_equat         = 0.0_dp       !< latitude of rotated equator of COSMO-DE grid [rad]
     80    REAL(dp) ::  phi_n             = 0.0_dp       !< latitude of rotated pole of COSMO-DE grid [rad]
     81    REAL(dp) ::  lambda_n          = 0.0_dp       !< longitude of rotaded pole of COSMO-DE grid [rad]
     82    REAL(dp) ::  phi_c             = 0.0_dp       !< rotated-grid latitude of the center of the PALM domain [rad]
     83    REAL(dp) ::  lambda_c          = 0.0_dp       !< rotated-grid longitude of the centre of the PALM domain [rad]
     84    REAL(dp) ::  phi_cn            = 0.0_dp       !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
     85    REAL(dp) ::  lambda_cn         = 0.0_dp       !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
     86    REAL(dp) ::  gam               = 0.0_dp       !< angle for working around phirot2phi/rlarot2rla bug
     87    REAL(dp) ::  dx                = 0.0_dp       !< PALM-4U grid spacing in x direction [m]
     88    REAL(dp) ::  dy                = 0.0_dp       !< PALM-4U grid spacing in y direction [m]
     89    REAL(dp) ::  dz(10)            = -1.0_dp      !< PALM-4U grid spacing in z direction [m]
     90    REAL(dp) ::  dz_max            = 1000.0_dp    !< maximum vertical grid spacing [m]
     91    REAL(dp) ::  dz_stretch_factor = 1.08_dp      !< factor for vertical grid stretching [m]
     92    REAL(dp) ::  dz_stretch_level  = -9999999.9_dp!< height above which the vertical grid will be stretched [m]
     93    REAL(dp) ::  dz_stretch_level_start(9) = -9999999.9_dp !< namelist parameter
     94    REAL(dp) ::  dz_stretch_level_end(9) = 9999999.9_dp !< namelist parameter
     95    REAL(dp) ::  dz_stretch_factor_array(9) = 1.08_dp !< namelist parameter
     96    REAL(dp) ::  dxi               = 0.0_dp       !< inverse PALM-4U grid spacing in x direction [m^-1]
     97    REAL(dp) ::  dyi               = 0.0_dp       !< inverse PALM-4U grid spacing in y direction [m^-1]
     98    REAL(dp) ::  dzi               = 0.0_dp       !< inverse PALM-4U grid spacing in z direction [m^-1]
     99    REAL(dp) ::  lx                = 0.0_dp       !< PALM-4U domain size in x direction [m]
     100    REAL(dp) ::  ly                = 0.0_dp       !< PALM-4U domain size in y direction [m]
     101    REAL(dp) ::  ug                = 0.0_dp       !< geostrophic wind in x direction [m/s]
     102    REAL(dp) ::  vg                = 0.0_dp       !< geostrophic wind in y direction [m/s]
     103    REAL(dp) ::  p0                = 0.0_dp       !< PALM-4U surface pressure, at z0 [Pa]
     104    REAL(dp) ::  x0                = 0.0_dp       !< x coordinate of PALM-4U Earth tangent [m]
     105    REAL(dp) ::  y0                = 0.0_dp       !< y coordinate of PALM-4U Earth tangent [m]
     106    REAL(dp) ::  z0                = 0.0_dp       !< Elevation of the PALM-4U domain above sea level [m]
     107    REAL(dp) ::  z_top             = 0.0_dp       !< height of the scalar top boundary [m]
     108    REAL(dp) ::  zw_top            = 0.0_dp       !< height of the vertical velocity top boundary [m]
     109    REAL(dp) ::  lonmin            = 0.0_dp       !< Minimunm longitude of COSMO-DE's rotated-pole grid
     110    REAL(dp) ::  lonmax            = 0.0_dp       !< Maximum longitude of COSMO-DE's rotated-pole grid
     111    REAL(dp) ::  latmin            = 0.0_dp       !< Minimunm latitude of COSMO-DE's rotated-pole grid
     112    REAL(dp) ::  latmax            = 0.0_dp       !< Maximum latitude of COSMO-DE's rotated-pole grid
     113    REAL(dp) ::  latitude          = 0.0_dp       !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
     114    REAL(dp) ::  longitude         = 0.0_dp       !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
     115    REAL(dp) ::  origin_lat        = 0.0_dp       !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg]
     116    REAL(dp) ::  origin_lon        = 0.0_dp       !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg]
     117    REAL(dp) ::  rotation_angle    = 0.0_dp       !< clockwise angle the PALM-4U north is rotated away from geographical north [deg]
     118    REAL(dp) ::  end_time          = 0.0_dp       !< PALM-4U simulation time [s]
    101119
    102120    REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hhl             !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
     
    106124    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlon            !< longitudes of COSMO-DE's rotated-pole grid
    107125    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlat            !< latitudes of COSMO-DE's rotated-pole grid
    108     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  time
     126    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  time            !< output times
     127    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  x               !< base palm grid x coordinate vector pointed to by grid_definitions
     128    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  xu              !< base palm grid xu coordinate vector pointed to by grid_definitions
     129    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  y               !< base palm grid y coordinate vector pointed to by grid_definitions
     130    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  yv              !< base palm grid yv coordinate vector pointed to by grid_definitions
     131    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z_column        !< base palm grid z coordinate vector including the top boundary coordinate (entire column)
     132    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw_column       !< base palm grid zw coordinate vector including the top boundary coordinate (entire column)
     133    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z               !< base palm grid z coordinate vector pointed to by grid_definitions
     134    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw              !< base palm grid zw coordinate vector pointed to by grid_definitions
    109135
    110136    INTEGER(hp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  soiltyp      !< COSMO-DE soil type map
     137    INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
     138    INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
    111139    INTEGER ::  i     !< indexing variable
    112     INTEGER ::  imin, imax, jmin, jmax !< index ranges for profile averaging
     140    INTEGER ::  average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging
    113141    INTEGER ::  k     !< indexing variable
    114142    INTEGER ::  nt    !< number of output time steps
     
    130158    LOGICAL ::  init_variables_required
    131159    LOGICAL ::  boundary_variables_required
     160    LOGICAL ::  ls_forcing_variables_required
     161    LOGICAL ::  profile_grids_required
    132162
    133163    INTEGER ::  n_invar = 0 !< number of variables in the input variable table
     
    193223    TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
    194224 
    195     NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude
     225    NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude,             &
     226                      dz_max, dz_stretch_factor, dz_stretch_level,             & !< old grid stretching parameters
     227                      dz_stretch_level_start, dz_stretch_level_end               !< new grid stretching parameters
    196228    NAMELIST /d3par/  end_time
    197229   
    198230    CHARACTER(LEN=LNAME) ::  nc_source_text     = ''  !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...'
    199     CHARACTER(LEN=DATE)  ::  start_date         = ''  !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation
    200     CHARACTER(LEN=PATH)  ::  hhl_file           = ''  !< Path to the file containing the COSMO-DE HHL variable (height of half layers, i.e. vertical cell faces)
    201     CHARACTER(LEN=PATH)  ::  namelist_file      = ''  !< Path to the PALM-4U namelist file
    202     CHARACTER(LEN=PATH)  ::  static_driver_file = ''  !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
    203     CHARACTER(LEN=PATH)  ::  soiltyp_file       = ''  !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
    204     CHARACTER(LEN=PATH)  ::  input_path         = ''  !< Path to the input data file directory
    205231    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  flow_files
    206232    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_moisture_files
     
    212238    CHARACTER(LEN=SNAME) ::  radiation_suffix     !< Suffix of radiation input files, e.g. 'radiation'
    213239    CHARACTER(LEN=SNAME) ::  soilmoisture_suffix  !< Suffix of input files for soil moisture spin-up, e.g. 'soilmoisture'
    214     CHARACTER(LEN=SNAME) ::  mode                 !< INIFOR's initialization mode, 'profile' or 'volume'
    215240                         
    216241    TYPE(nc_file) ::  output_file
     242
     243    TYPE(inifor_config) ::  cfg !< Container of INIFOR configuration
    217244
    218245 CONTAINS
     
    224251! Section 1: Define default parameters
    225252!------------------------------------------------------------------------------
    226        start_date = '2013072100'
     253       cfg % start_date = '2013072100'
    227254       end_hour = 2
    228255       start_hour_soil = -2
     
    249276       origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
    250277       origin_lon = 13.082744_dp * TO_RADIANS
    251        z0 = 35.0_dp
     278       cfg % z0 = 35.0_dp
    252279
    253280       ! Default atmospheric parameters
    254281       ug = 0.0_dp
    255282       vg = 0.0_dp
    256        p0 = P_SL
     283       cfg % p0 = P_SL
    257284
    258285       ! Parameters for file names
     
    262289       start_hour_soilmoisture = start_hour_flow - 2
    263290       end_hour_soilmoisture = start_hour_flow
    264        step_hour = 1
     291       step_hour = FORCING_STEP
     292
    265293       input_prefix = 'laf'  ! 'laf' for COSMO-DE analyses
     294       cfg % flow_prefix = input_prefix
     295       cfg % soil_prefix = input_prefix
     296       cfg % radiation_prefix = input_prefix
     297       cfg % soilmoisture_prefix  = input_prefix
     298
    266299       flow_suffix = '-flow'
    267300       soil_suffix = '-soil'
     
    273306!------------------------------------------------------------------------------
    274307
    275        ! Set default paths
    276        input_path         = './'
    277        hhl_file           = ''
    278        soiltyp_file       = ''
    279        namelist_file      = './namelist'
    280        output_file % name = './palm-4u-input.nc'
     308       ! Set default paths and modes
     309       cfg % input_path         = './'
     310       cfg % hhl_file           = ''
     311       cfg % soiltyp_file       = ''
     312       cfg % namelist_file      = './namelist'
     313       cfg % static_driver_file = ''
     314       cfg % output_file = './palm-4u-input.nc'
     315       cfg % ic_mode = 'volume'
     316       cfg % bc_mode = 'real'
    281317
    282318       ! Update default file names and domain centre
    283        CALL parse_command_line_arguments( start_date, hhl_file, soiltyp_file,  &
    284                static_driver_file, input_path, output_file % name,             &
    285                namelist_file, ug, vg, p0, z0, mode )
    286 
    287        init_variables_required     = .TRUE.
    288        boundary_variables_required = (TRIM(mode) .NE. 'profile')
    289 
    290        CALL normalize_path(input_path)
    291        IF (TRIM(hhl_file) == '')  hhl_file = TRIM(input_path) // 'hhl.nc'
    292        IF (TRIM(soiltyp_file) == '')  soiltyp_file = TRIM(input_path) // 'soil.nc'
    293 
    294        CALL report('setup_parameters', "       data path: " // TRIM(input_path))
    295        CALL report('setup_parameters', "        hhl file: " // TRIM(hhl_file))
    296        CALL report('setup_parameters', "    soiltyp file: " // TRIM(soiltyp_file))
    297        CALL report('setup_parameters', "   namelist file: " // TRIM(namelist_file))
    298        CALL report('setup_parameters', "output data file: " // TRIM(output_file % name))
     319       CALL parse_command_line_arguments( cfg )
     320
     321       output_file % name = cfg % output_file
     322       z0 = cfg % z0
     323       p0 = cfg % p0
     324
     325       init_variables_required = .TRUE.
     326       boundary_variables_required = TRIM( cfg % bc_mode ) == 'real'
     327       ls_forcing_variables_required = TRIM( cfg % bc_mode ) == 'ideal'
     328
     329       IF ( ls_forcing_variables_required )  THEN
     330          message = "Averaging of large-scale forcing profiles " //            &
     331                    "has not been implemented, yet."
     332          CALL abort('setup_parameters', message)
     333       END IF
     334
     335       CALL normalize_path(cfg % input_path)
     336       IF (TRIM(cfg % hhl_file) == '')  cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
     337       IF (TRIM(cfg % soiltyp_file) == '')  cfg % soiltyp_file = TRIM(cfg % input_path) // 'soil.nc'
     338
     339       CALL validate_config( cfg )
     340
     341       CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode))
     342       CALL report('setup_parameters', "       forcing mode: " // TRIM(cfg % bc_mode))
     343       CALL report('setup_parameters', "          data path: " // TRIM(cfg % input_path))
     344       CALL report('setup_parameters', "           hhl file: " // TRIM(cfg % hhl_file))
     345       CALL report('setup_parameters', "       soiltyp file: " // TRIM(cfg % soiltyp_file))
     346       CALL report('setup_parameters', "      namelist file: " // TRIM(cfg % namelist_file))
     347       CALL report('setup_parameters', "   output data file: " // TRIM(output_file % name))
    299348
    300349 CALL run_control('time', 'init')
    301350       ! Read in namelist parameters
    302        OPEN(10, FILE=namelist_file)
     351       OPEN(10, FILE=cfg % namelist_file)
    303352       READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz
    304353       READ(10, NML=d3par)  ! end_time
     
    306355 CALL run_control('time', 'read')
    307356
    308        end_hour = CEILING(end_time / FORCING_FREQ)
     357       end_hour = CEILING( end_time / 3600.0 * step_hour )
    309358
    310359       ! Generate input file lists
    311        CALL input_file_list(start_date, start_hour_flow, end_hour, step_hour,  &
    312           input_path, input_prefix, flow_suffix, flow_files)
    313        CALL input_file_list(start_date, start_hour_soil, end_hour, step_hour,  &
    314           input_path, input_prefix, soil_suffix, soil_files)
    315        CALL input_file_list(start_date, start_hour_radiation, end_hour, step_hour, &
    316           input_path, input_prefix, radiation_suffix, radiation_files)
    317        CALL input_file_list(start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
    318           input_path, input_prefix, soilmoisture_suffix, soil_moisture_files)
     360       CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour,  &
     361          cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files)
     362       CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour,  &
     363          cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files)
     364       CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, &
     365          cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files)
     366       CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
     367          cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
    319368
    320369!
     
    322371! Section 3: Check for consistency
    323372!------------------------------------------------------------------------------
    324        IF (dx*dy*dz .EQ. 0.0_dp)  THEN
    325           message = "Grid cells have zero volume. Grid spacings are probably"//&
    326              " set incorrectly in namelist file '" // TRIM(namelist_file) // "'."
    327           CALL abort('setup_parameters', message)
    328        END IF
     373
    329374!
    330375!------------------------------------------------------------------------------
     
    339384       ! Read COSMO-DE soil type map
    340385       cosmo_var % name = 'SOILTYP'
    341        soiltyp = NINT(get_netcdf_variable(soiltyp_file, cosmo_var), hp)
     386       CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp)
    342387
    343388       message = 'Reading PALM-4U origin from'
    344        IF (TRIM(static_driver_file) .NE. '')  THEN
    345 
    346           origin_lon = get_netcdf_attribute(static_driver_file, 'origin_lon')
    347           origin_lat = get_netcdf_attribute(static_driver_file, 'origin_lat')
     389       IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
     390
     391          origin_lon = get_netcdf_attribute(cfg % static_driver_file, 'origin_lon')
     392          origin_lat = get_netcdf_attribute(cfg % static_driver_file, 'origin_lat')
    348393
    349394          message = TRIM(message) // " static driver file '"                   &
    350                                   // TRIM(static_driver_file) // "'"
     395                                  // TRIM(cfg % static_driver_file) // "'"
    351396
    352397
     
    357402
    358403          message = TRIM(message) // " namlist file '"                         &
    359                                   // TRIM(namelist_file) // "'"
     404                                  // TRIM(cfg % namelist_file) // "'"
    360405
    361406       END IF
     
    368413       ! Read in COSMO-DE heights of half layers (vertical cell faces)
    369414       cosmo_var % name = 'HHL'
    370        hhl = get_netcdf_variable(hhl_file, cosmo_var)
     415       CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl)
    371416 CALL run_control('time', 'read')
    372417
     
    392437       lx = (nx+1) * dx
    393438       ly = (ny+1) * dy
    394        lz = (nz+1) * dz
    395439       
    396440       ! PALM-4U point of Earth tangency
     
    399443
    400444       ! time vector
    401        nt = CEILING(end_time / FORCING_FREQ) + 1
     445       nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1
    402446       ALLOCATE( time(nt) )
    403447       CALL linspace(0.0_dp, 3600.0_dp * (nt-1), time)
     
    463507    SUBROUTINE setup_grids() ! setup_grids(inifor_settings(with nx, ny, nz,...))
    464508       CHARACTER ::  interp_mode
    465        
     509
    466510!------------------------------------------------------------------------------
    467 ! Section 1: Define model and initialization grids
     511! Section 0: Define base PALM-4U grid coordinate vectors
     512!------------------------------------------------------------------------------
     513       ! palm x y z, we allocate the column to nz+1 in order to include the top
     514       ! scalar boundary. The interpolation grids will be associated with
     515       ! a shorter column that omits the top element.
     516       
     517       ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
     518       CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x)
     519       CALL linspace(0.5_dp * dy, ly - 0.5_dp * dy, y)
     520       CALL stretched_z(z_column, dz, dz_max=dz_max,                           &
     521                        dz_stretch_factor=dz_stretch_factor,                   &
     522                        dz_stretch_level=dz_stretch_level,                     &
     523                        dz_stretch_level_start=dz_stretch_level_start,         &
     524                        dz_stretch_level_end=dz_stretch_level_end,             &
     525                        dz_stretch_factor_array=dz_stretch_factor_array)
     526       z(1:nz) = z_column(1:nz)
     527       z_top = z_column(nz+1)
     528
     529       ! palm xu yv zw, compared to the scalar grid, velocity coordinates
     530       ! contain one element less.
     531       ALLOCATE( xu(1:nx),  yv(1:ny), zw(1:nz-1), zw_column(1:nz))
     532       CALL linspace(dx, lx - dx, xu)
     533       CALL linspace(dy, ly - dy, yv)
     534       CALL midpoints(z_column, zw_column)
     535       zw(1:nz-1) = zw_column(1:nz-1)
     536       zw_top     = zw_column(nz)
     537
     538
     539!------------------------------------------------------------------------------
     540! Section 1: Define initialization and boundary grids
    468541!------------------------------------------------------------------------------
    469542       CALL init_grid_definition('palm', grid=palm_grid,                       &
    470543               xmin=0.0_dp, xmax=lx,                                           &
    471544               ymin=0.0_dp, ymax=ly,                                           &
    472                zmin=0.0_dp, zmax=lz, x0=x0, y0=y0, z0=z0,                      &
    473                nx=nx, ny=ny, nz=nz, mode=mode)
     545               x0=x0, y0=y0, z0=z0,                      &
     546               nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode)
    474547
    475548       ! Subtracting 1 because arrays will be allocated with nlon + 1 elements.
     
    477550               xmin=lonmin, xmax=lonmax,                                       &
    478551               ymin=latmin, ymax=latmax,                                       &
    479                zmin=0.0_dp, zmax=51.0_dp, x0=x0, y0=y0, z0=0.0_dp,             &
     552               x0=x0, y0=y0, z0=0.0_dp,             &
    480553               nx=nlon-1, ny=nlat-1, nz=nlev-1)
    481554
     
    487560               xmin=0.0_dp, xmax=lx,                                           &
    488561               ymin=0.0_dp, ymax=ly,                                           &
    489                zmin=0.0_dp, zmax=lz, x0=x0, y0=y0, z0=z0,                      &
     562               x0=x0, y0=y0, z0=z0,                      &
    490563               nx=nx, ny=ny, nz=nlev-2)
    491564
     
    493566               xmin = dx, xmax = lx - dx,                                      &
    494567               ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    495                zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    496568               x0=x0, y0=y0, z0 = z0,                                          &
    497569               nx = nx-1, ny = ny, nz = nz,                                    &
    498                dx = dx, dy = dy, dz = dz, mode=mode)
     570               z=z, ic_mode=cfg % ic_mode)
    499571
    500572       CALL init_grid_definition('boundary', grid=v_initial_grid,              &
    501573               xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    502574               ymin = dy, ymax = ly - dy,                                      &
    503                zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    504575               x0=x0, y0=y0, z0 = z0,                                          &
    505576               nx = nx, ny = ny-1, nz = nz,                                    &
    506                dx = dx, dy = dy, dz = dz, mode=mode)
     577               z=z, ic_mode=cfg % ic_mode)
    507578
    508579       CALL init_grid_definition('boundary', grid=w_initial_grid,              &
    509580               xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    510581               ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    511                zmin = dz, zmax = lz - dz,                                      &
    512582               x0=x0, y0=y0, z0 = z0,                                          &
    513583               nx = nx, ny = ny, nz = nz-1,                                    &
    514                dx = dx, dy = dy, dz = dz, mode=mode)
     584               z=zw, ic_mode=cfg % ic_mode)
    515585
    516586       CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate,      &
    517587               xmin = dx, xmax = lx - dx,                                      &
    518588               ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    519                zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    520589               x0=x0, y0=y0, z0 = z0,                                          &
    521                nx = nx-1, ny = ny, nz = nlev - 2,                              &
    522                dx = dx, dy = dy, dz = dz)
     590               nx = nx-1, ny = ny, nz = nlev - 2)
    523591
    524592       CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate,      &
    525593               xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    526594               ymin = dy, ymax = ly - dy,                                      &
    527                zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    528595               x0=x0, y0=y0, z0 = z0,                                          &
    529                nx = nx, ny = ny-1, nz = nlev - 2,                              &
    530                dx = dx, dy = dy, dz = dz)
     596               nx = nx, ny = ny-1, nz = nlev - 2)
    531597
    532598       CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate,      &
    533599               xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    534600               ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    535                zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    536601               x0=x0, y0=y0, z0 = z0,                                          &
    537                nx = nx, ny = ny, nz = nlev - 2,                                &
    538                dx = dx, dy = dy, dz = dz)
     602               nx = nx, ny = ny, nz = nlev - 2)
    539603
    540604      IF (boundary_variables_required)  THEN
     
    546610                  xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    547611                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    548                   zmin =  0.5_dp * dz, zmax = lz - 0.5_dp * dz,                   &
    549                   x0=x0, y0=y0, z0 = z0,                                          &
    550                   nx = 0, ny = ny, nz = nz,                                       &
    551                   dx = dx, dy = dy, dz = dz)
     612                  x0=x0, y0=y0, z0 = z0,                                          &
     613                  nx = 0, ny = ny, nz = nz, z=z)
    552614
    553615          CALL init_grid_definition('boundary', grid=scalars_west_grid,           &
    554616                  xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    555617                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    556                   zmin =  0.5_dp * dz, zmax = lz - 0.5_dp * dz,                   &
    557                   x0=x0, y0=y0, z0 = z0,                                          &
    558                   nx = 0, ny = ny, nz = nz,                                       &
    559                   dx = dx, dy = dy, dz = dz)
     618                  x0=x0, y0=y0, z0 = z0,                                          &
     619                  nx = 0, ny = ny, nz = nz, z=z)
    560620
    561621          CALL init_grid_definition('boundary', grid=scalars_north_grid,          &
    562622                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    563623                  ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    564                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    565                   x0=x0, y0=y0, z0 = z0,                                          &
    566                   nx = nx, ny = 0, nz = nz,                                       &
    567                   dx = dx, dy = dy, dz = dz)
     624                  x0=x0, y0=y0, z0 = z0,                                          &
     625                  nx = nx, ny = 0, nz = nz, z=z)
    568626
    569627          CALL init_grid_definition('boundary', grid=scalars_south_grid,          &
    570628                  xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    571629                  ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    572                   zmin =  0.5_dp * dz, zmax = lz - 0.5_dp * dz,                   &
    573                   x0=x0, y0=y0, z0 = z0,                                          &
    574                   nx = nx, ny = 0, nz = nz,                                       &
    575                   dx = dx, dy = dy, dz = dz)
     630                  x0=x0, y0=y0, z0 = z0,                                          &
     631                  nx = nx, ny = 0, nz = nz, z=z)
    576632
    577633          CALL init_grid_definition('boundary', grid=scalars_top_grid,            &
    578634                  xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    579635                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    580                   zmin =  lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz,              &
    581                   x0=x0, y0=y0, z0 = z0,                                          &
    582                   nx = nx, ny = ny, nz = 0,                                       &
    583                   dx = dx, dy = dy, dz = dz)
     636                  x0=x0, y0=y0, z0 = z0,                                          &
     637                  nx = nx, ny = ny, nz = 1, z=(/z_top/))
    584638
    585639          CALL init_grid_definition('boundary', grid=u_east_grid,                 &
    586640                  xmin = lx, xmax = lx,                                           &
    587641                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    588                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    589                   x0=x0, y0=y0, z0 = z0,                                          &
    590                   nx = 0, ny = ny, nz = nz,                                       &
    591                   dx = dx, dy = dy, dz = dz)
     642                  x0=x0, y0=y0, z0 = z0,                                          &
     643                  nx = 0, ny = ny, nz = nz, z=z)
    592644
    593645          CALL init_grid_definition('boundary', grid=u_west_grid,                 &
    594646                  xmin = 0.0_dp, xmax = 0.0_dp,                                   &
    595647                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    596                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    597                   x0=x0, y0=y0, z0 = z0,                                          &
    598                   nx = 0, ny = ny, nz = nz,                                       &
    599                   dx = dx, dy = dy, dz = dz)
     648                  x0=x0, y0=y0, z0 = z0,                                          &
     649                  nx = 0, ny = ny, nz = nz, z=z)
    600650
    601651          CALL init_grid_definition('boundary', grid=u_north_grid,                &
    602652                  xmin = dx, xmax = lx - dx,                                      &
    603653                  ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    604                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    605                   x0=x0, y0=y0, z0 = z0,                                          &
    606                   nx = nx-1, ny = 0, nz = nz,                                     &
    607                   dx = dx, dy = dy, dz = dz)
    608 
     654                  x0=x0, y0=y0, z0 = z0,                                          &
     655                  nx = nx-1, ny = 0, nz = nz, z=z)
     656   
    609657          CALL init_grid_definition('boundary', grid=u_south_grid,                &
    610658                  xmin = dx, xmax = lx - dx,                                      &
    611659                  ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    612                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    613                   x0=x0, y0=y0, z0 = z0,                                          &
    614                   nx = nx-1, ny = 0, nz = nz,                                     &
    615                   dx = dx, dy = dy, dz = dz)
     660                  x0=x0, y0=y0, z0 = z0,                                          &
     661                  nx = nx-1, ny = 0, nz = nz, z=z)
    616662
    617663          CALL init_grid_definition('boundary', grid=u_top_grid,                  &
    618664                  xmin = dx, xmax = lx - dx,                                      &
    619665                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    620                   zmin = lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz,               &
    621                   x0=x0, y0=y0, z0 = z0,                                          &
    622                   nx = nx-1, ny = ny, nz = 0,                                     &
    623                   dx = dx, dy = dy, dz = dz)
     666                  x0=x0, y0=y0, z0 = z0,                                          &
     667                  nx = nx-1, ny = ny, nz = 1, z=(/z_top/))
    624668
    625669          CALL init_grid_definition('boundary', grid=v_east_grid,                 &
    626670                  xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    627671                  ymin = dy, ymax = ly - dy,                                      &
    628                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    629                   x0=x0, y0=y0, z0 = z0,                                          &
    630                   nx = 0, ny = ny-1, nz = nz,                                     &
    631                   dx = dx, dy = dy, dz = dz)
     672                  x0=x0, y0=y0, z0 = z0,                                          &
     673                  nx = 0, ny = ny-1, nz = nz, z=z)
    632674
    633675          CALL init_grid_definition('boundary', grid=v_west_grid,                 &
    634676                  xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    635677                  ymin = dy, ymax = ly - dy,                                      &
    636                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    637                   x0=x0, y0=y0, z0 = z0,                                          &
    638                   nx = 0, ny = ny-1, nz = nz,                                     &
    639                   dx = dx, dy = dy, dz = dz)
     678                  x0=x0, y0=y0, z0 = z0,                                          &
     679                  nx = 0, ny = ny-1, nz = nz, z=z)
    640680
    641681          CALL init_grid_definition('boundary', grid=v_north_grid,                &
    642682                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    643683                  ymin = ly, ymax = ly,                                           &
    644                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    645                   x0=x0, y0=y0, z0 = z0,                                          &
    646                   nx = nx, ny = 0, nz = nz,                                       &
    647                   dx = dx, dy = dy, dz = dz)
     684                  x0=x0, y0=y0, z0 = z0,                                          &
     685                  nx = nx, ny = 0, nz = nz, z=z)
    648686
    649687          CALL init_grid_definition('boundary', grid=v_south_grid,                &
    650688                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    651689                  ymin = 0.0_dp, ymax = 0.0_dp,                                   &
    652                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
    653                   x0=x0, y0=y0, z0 = z0,                                          &
    654                   nx = nx, ny = 0, nz = nz,                                       &
    655                   dx = dx, dy = dy, dz = dz)
     690                  x0=x0, y0=y0, z0 = z0,                                          &
     691                  nx = nx, ny = 0, nz = nz, z=z)
    656692
    657693          CALL init_grid_definition('boundary', grid=v_top_grid,                  &
    658694                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    659695                  ymin = dy, ymax = ly - dy,                                      &
    660                   zmin = lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz,               &
    661                   x0=x0, y0=y0, z0 = z0,                                          &
    662                   nx = nx, ny = ny-1, nz = 0,                                     &
    663                   dx = dx, dy = dy, dz = dz)
     696                  x0=x0, y0=y0, z0 = z0,                                          &
     697                  nx = nx, ny = ny-1, nz = 1, z=(/z_top/))
    664698
    665699          CALL init_grid_definition('boundary', grid=w_east_grid,                 &
    666700                  xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    667701                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    668                   zmin =  dz, zmax = lz - dz,                                     &
    669                   x0=x0, y0=y0, z0 = z0,                                          &
    670                   nx = 0, ny = ny, nz = nz - 1,                                   &
    671                   dx = dx, dy = dy, dz = dz)
     702                  x0=x0, y0=y0, z0 = z0,                                          &
     703                  nx = 0, ny = ny, nz = nz - 1, z=zw)
    672704
    673705          CALL init_grid_definition('boundary', grid=w_west_grid,                 &
    674706                  xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    675707                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    676                   zmin = dz, zmax = lz - dz,                                      &
    677                   x0=x0, y0=y0, z0 = z0,                                          &
    678                   nx = 0, ny = ny, nz = nz - 1,                                   &
    679                   dx = dx, dy = dy, dz = dz)
     708                  x0=x0, y0=y0, z0 = z0,                                          &
     709                  nx = 0, ny = ny, nz = nz - 1, z=zw)
    680710
    681711          CALL init_grid_definition('boundary', grid=w_north_grid,                &
    682712                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    683713                  ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    684                   zmin = dz, zmax = lz - dz,                                      &
    685                   x0=x0, y0=y0, z0 = z0,                                          &
    686                   nx = nx, ny = 0, nz = nz - 1,                                   &
    687                   dx = dx, dy = dy, dz = dz)
     714                  x0=x0, y0=y0, z0 = z0,                                          &
     715                  nx = nx, ny = 0, nz = nz - 1, z=zw)
    688716
    689717          CALL init_grid_definition('boundary', grid=w_south_grid,                &
    690718                  xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    691719                  ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    692                   zmin = dz, zmax = lz - dz,                                      &
    693                   x0=x0, y0=y0, z0 = z0,                                          &
    694                   nx = nx, ny = 0, nz = nz - 1,                                   &
    695                   dx = dx, dy = dy, dz = dz)
     720                  x0=x0, y0=y0, z0 = z0,                                          &
     721                  nx = nx, ny = 0, nz = nz - 1, z=zw)
    696722
    697723          CALL init_grid_definition('boundary', grid=w_top_grid,                  &
    698724                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    699725                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    700                   zmin = lz, zmax = lz,                                           &
    701                   x0=x0, y0=y0, z0 = z0,                                          &
    702                   nx = nx, ny = ny, nz = 0,                                       &
    703                   dx = dx, dy = dy, dz = dz)
     726                  x0=x0, y0=y0, z0 = z0,                                          &
     727                  nx = nx, ny = ny, nz = 1, z=(/zw_top/))
    704728
    705729          CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate,   &
    706730                  xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    707731                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    708                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    709                   x0=x0, y0=y0, z0 = z0,                                          &
    710                   nx = 0, ny = ny, nz = nlev - 2,                                 &
    711                   dx = dx, dy = dy, dz = dz)
     732                  x0=x0, y0=y0, z0 = z0,                                          &
     733                  nx = 0, ny = ny, nz = nlev - 2)
    712734
    713735          CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate,   &
    714736                  xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    715737                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    716                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    717                   x0=x0, y0=y0, z0 = z0,                                          &
    718                   nx = 0, ny = ny, nz = nlev - 2,                                 &
    719                   dx = dx, dy = dy, dz = dz)
     738                  x0=x0, y0=y0, z0 = z0,                                          &
     739                  nx = 0, ny = ny, nz = nlev - 2)
    720740
    721741          CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate,  &
    722742                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    723743                  ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    724                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    725                   x0=x0, y0=y0, z0 = z0,                                          &
    726                   nx = nx, ny = 0, nz = nlev - 2,                                 &
    727                   dx = dx, dy = dy, dz = dz)
     744                  x0=x0, y0=y0, z0 = z0,                                          &
     745                  nx = nx, ny = 0, nz = nlev - 2)
    728746
    729747          CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate,  &
    730748                  xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    731749                  ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    732                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    733                   x0=x0, y0=y0, z0 = z0,                                          &
    734                   nx = nx, ny = 0, nz = nlev - 2,                                 &
    735                   dx = dx, dy = dy, dz = dz)
     750                  x0=x0, y0=y0, z0 = z0,                                          &
     751                  nx = nx, ny = 0, nz = nlev - 2)
    736752
    737753          CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate,    &
    738754                  xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    739755                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    740                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    741                   x0=x0, y0=y0, z0 = z0,                                          &
    742                   nx = nx, ny = ny, nz = nlev - 2,                                &
    743                   dx = dx, dy = dy, dz = dz)
     756                  x0=x0, y0=y0, z0 = z0,                                          &
     757                  nx = nx, ny = ny, nz = nlev - 2)
    744758
    745759          CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate,         &
    746760                  xmin = lx, xmax = lx,                                           &
    747761                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    748                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    749                   x0=x0, y0=y0, z0 = z0,                                          &
    750                   nx = 0, ny = ny, nz = nlev - 2,                                 &
    751                   dx = dx, dy = dy, dz = dz)
     762                  x0=x0, y0=y0, z0 = z0,                                          &
     763                  nx = 0, ny = ny, nz = nlev - 2)
    752764
    753765          CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate,         &
    754766                  xmin = 0.0_dp, xmax = 0.0_dp,                                   &
    755767                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    756                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    757                   x0=x0, y0=y0, z0 = z0,                                          &
    758                   nx = 0, ny = ny, nz = nlev - 2,                                 &
    759                   dx = dx, dy = dy, dz = dz)
     768                  x0=x0, y0=y0, z0 = z0,                                          &
     769                  nx = 0, ny = ny, nz = nlev - 2)
    760770
    761771          CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate,        &
    762772                  xmin = dx, xmax = lx - dx,                                      &
    763773                  ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    764                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    765                   x0=x0, y0=y0, z0 = z0,                                          &
    766                   nx = nx-1, ny = 0, nz = nlev - 2,                               &
    767                   dx = dx, dy = dy, dz = dz)
     774                  x0=x0, y0=y0, z0 = z0,                                          &
     775                  nx = nx-1, ny = 0, nz = nlev - 2)
    768776
    769777          CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate,        &
    770778                  xmin = dx, xmax = lx - dx,                                      &
    771779                  ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    772                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    773                   x0=x0, y0=y0, z0 = z0,                                          &
    774                   nx = nx-1, ny = 0, nz = nlev - 2,                               &
    775                   dx = dx, dy = dy, dz = dz)
     780                  x0=x0, y0=y0, z0 = z0,                                          &
     781                  nx = nx-1, ny = 0, nz = nlev - 2)
    776782
    777783          CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate,          &
    778784                  xmin = dx, xmax = lx - dx,                                      &
    779785                  ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    780                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    781                   x0=x0, y0=y0, z0 = z0,                                          &
    782                   nx = nx-1, ny = ny, nz = nlev - 2,                              &
    783                   dx = dx, dy = dy, dz = dz)
     786                  x0=x0, y0=y0, z0 = z0,                                          &
     787                  nx = nx-1, ny = ny, nz = nlev - 2)
    784788
    785789          CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate,         &
    786790                  xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    787791                  ymin = dy, ymax = ly - dy,                                      &
    788                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    789                   x0=x0, y0=y0, z0 = z0,                                          &
    790                   nx = 0, ny = ny-1, nz = nlev - 2,                               &
    791                   dx = dx, dy = dy, dz = dz)
     792                  x0=x0, y0=y0, z0 = z0,                                          &
     793                  nx = 0, ny = ny-1, nz = nlev - 2)
    792794
    793795          CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate,         &
    794796                  xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    795797                  ymin = dy, ymax = ly - dy,                                      &
    796                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    797                   x0=x0, y0=y0, z0 = z0,                                          &
    798                   nx = 0, ny = ny-1, nz = nlev - 2,                               &
    799                   dx = dx, dy = dy, dz = dz)
     798                  x0=x0, y0=y0, z0 = z0,                                          &
     799                  nx = 0, ny = ny-1, nz = nlev - 2)
    800800
    801801          CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate,        &
    802802                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    803803                  ymin = ly, ymax = ly,                                           &
    804                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    805                   x0=x0, y0=y0, z0 = z0,                                          &
    806                   nx = nx, ny = 0, nz = nlev - 2,                                 &
    807                   dx = dx, dy = dy, dz = dz)
     804                  x0=x0, y0=y0, z0 = z0,                                          &
     805                  nx = nx, ny = 0, nz = nlev - 2)
    808806
    809807          CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate,        &
    810808                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    811809                  ymin = 0.0_dp, ymax = 0.0_dp,                                   &
    812                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    813                   x0=x0, y0=y0, z0 = z0,                                          &
    814                   nx = nx, ny = 0, nz = nlev - 2,                                 &
    815                   dx = dx, dy = dy, dz = dz)
     810                  x0=x0, y0=y0, z0 = z0,                                          &
     811                  nx = nx, ny = 0, nz = nlev - 2)
    816812
    817813          CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate,          &
    818814                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    819815                  ymin = dy, ymax = ly - dy,                                      &
    820                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    821                   x0=x0, y0=y0, z0 = z0,                                          &
    822                   nx = nx, ny = ny-1, nz = nlev - 2,                              &
    823                   dx = dx, dy = dy, dz = dz)
     816                  x0=x0, y0=y0, z0 = z0,                                          &
     817                  nx = nx, ny = ny-1, nz = nlev - 2)
    824818
    825819          CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate,         &
    826820                  xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    827821                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    828                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    829                   x0=x0, y0=y0, z0 = z0,                                          &
    830                   nx = 0, ny = ny, nz = nlev - 2,                                 &
    831                   dx = dx, dy = dy, dz = dz)
     822                  x0=x0, y0=y0, z0 = z0,                                          &
     823                  nx = 0, ny = ny, nz = nlev - 2)
    832824
    833825          CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
    834826                  xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    835827                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    836                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    837                   x0=x0, y0=y0, z0 = z0,                                          &
    838                   nx = 0, ny = ny, nz = nlev - 2,                                 &
    839                   dx = dx, dy = dy, dz = dz)
     828                  x0=x0, y0=y0, z0 = z0,                                          &
     829                  nx = 0, ny = ny, nz = nlev - 2)
    840830
    841831          CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
    842832                  xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    843833                  ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    844                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    845                   x0=x0, y0=y0, z0 = z0,                                          &
    846                   nx = nx, ny = 0, nz = nlev - 2,                                 &
    847                   dx = dx, dy = dy, dz = dz)
     834                  x0=x0, y0=y0, z0 = z0,                                          &
     835                  nx = nx, ny = 0, nz = nlev - 2)
    848836
    849837          CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
    850838                  xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    851839                  ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    852                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    853                   x0=x0, y0=y0, z0 = z0,                                          &
    854                   nx = nx, ny = 0, nz = nlev - 2,                                 &
    855                   dx = dx, dy = dy, dz = dz)
     840                  x0=x0, y0=y0, z0 = z0,                                          &
     841                  nx = nx, ny = 0, nz = nlev - 2)
    856842
    857843          CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
    858844                  xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    859845                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    860                   zmin =  0.0_dp, zmax = 0.0_dp,                                  &
    861                   x0=x0, y0=y0, z0 = z0,                                          &
    862                   nx = nx, ny = ny, nz = nlev - 2,                                &
    863                   dx = dx, dy = dy, dz = dz)
     846                  x0=x0, y0=y0, z0 = z0,                                          &
     847                  nx = nx, ny = ny, nz = nlev - 2)
    864848       END IF
    865849
     
    869853!------------------------------------------------------------------------------
    870854
    871        IF (TRIM(mode) == 'profile')  THEN
     855       profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR.              &
     856                                  TRIM(cfg % bc_mode) == 'ideal' )
     857
     858       IF (profile_grids_required)  THEN
    872859          CALL init_grid_definition('boundary', grid=scalar_profile_grid,      &
    873860                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    874861                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    875                   zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                 &
    876862                  x0=x0, y0=y0, z0 = z0,                                       &
    877                   nx = 0, ny = 0, nz = nz,                                     &
    878                   dx = dx, dy = dy, dz = dz)
     863                  nx = 0, ny = 0, nz = nz, z=z)
    879864
    880865          CALL init_grid_definition('boundary', grid=w_profile_grid,           &
    881866                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    882867                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    883                   zmin = dz, zmax = lz - dz,                                   &
    884868                  x0=x0, y0=y0, z0 = z0,                                       &
    885                   nx = 0, ny = 0, nz = nz - 1,                                 &
    886                   dx = dx, dy = dy, dz = dz)
     869                  nx = 0, ny = 0, nz = nz - 1, z=zw)
    887870
    888871          CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,&
    889872                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    890873                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    891                   zmin = 0.0_dp, zmax = 0.0_dp,                                &
    892874                  x0=x0, y0=y0, z0 = z0,                                       &
    893                   nx = 0, ny = 0, nz = nlev - 2,                               &
    894                   dx = dx, dy = dy, dz = dz)
     875                  nx = 0, ny = 0, nz = nlev - 2, z=z)
    895876
    896877          CALL init_grid_definition('boundary', grid=w_profile_intermediate,   &
    897878                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    898879                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    899                   zmin = 0.0_dp, zmax = 0.0_dp,                                &
    900880                  x0=x0, y0=y0, z0 = z0,                                       &
    901                   nx = 0, ny = 0, nz = nlev - 2,                               &
    902                   dx = dx, dy = dy, dz = dz)
     881                  nx = 0, ny = 0, nz = nlev - 2, z=zw)
    903882       END IF
    904883
     
    908887!------------------------------------------------------------------------------
    909888       interp_mode = 's'
    910        CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, mode=mode)
     889       CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    911890       IF (boundary_variables_required)  THEN
    912891          CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode)
     
    918897
    919898       interp_mode = 'u'
    920        CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, mode=mode)
     899       CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    921900       IF (boundary_variables_required)  THEN
    922901          CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode)
     
    928907
    929908       interp_mode = 'v'
    930        CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, mode=mode)
     909       CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    931910       IF (boundary_variables_required)  THEN
    932911          CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode)
     
    938917
    939918       interp_mode = 'w'
    940        CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, mode=mode)
     919       CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    941920       IF (boundary_variables_required)  THEN
    942921          CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode)
     
    947926       END IF
    948927
    949        IF (TRIM(mode) == 'profile')  THEN
    950            CALL setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax)
     928       IF (TRIM(cfg % ic_mode) == 'profile')  THEN
     929           CALL setup_averaging(cosmo_grid, palm_intermediate,                 &
     930                average_imin, average_imax, average_jmin, average_jmax)
    951931       END IF
    952932       
     
    955935
    956936
    957     SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, mode)
     937    SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode)
    958938
    959939       TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
    960940       TYPE(grid_definition), INTENT(INOUT), TARGET ::  palm_grid, palm_intermediate
    961941       CHARACTER, INTENT(IN)                        ::  kind
    962        CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  mode
     942       CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
    963943
    964944       TYPE(grid_definition), POINTER      ::  grid
     
    966946       REAL(dp), DIMENSION(:,:,:), POINTER ::  h
    967947
    968        LOGICAL :: setup_vertical = .TRUE.
    969 
    970        IF (PRESENT(mode))  THEN
    971           IF (TRIM(mode) == 'profile')  setup_vertical = .FALSE.
    972        ELSE
    973           setup_vertical = .TRUE.
     948       LOGICAL :: setup_vertical
     949
     950       setup_vertical = .TRUE.
     951       IF (PRESENT(ic_mode))  THEN
     952          IF (TRIM(ic_mode) == 'profile')  setup_vertical = .FALSE.
    974953       END IF
    975954
     
    1006985       CASE DEFAULT
    1007986
    1008           message = "Interpolation mode '" // mode // "' is not supported."
     987          message = "Interpolation quantity '" // kind // "' is not supported."
    1009988          CALL abort('setup_interpolation', message)
    1010989
     
    1014993
    1015994       CALL find_horizontal_neighbours(lat, lon,                               &
    1016           cosmo_grid % dxi, cosmo_grid % dyi, grid % clat,                     &
    1017           grid % clon, grid % ii, grid % jj)
     995          grid % clat, grid % clon, grid % ii, grid % jj)
    1018996
    1019997       CALL compute_horizontal_interp_weights(lat, lon,                        &
    1020           cosmo_grid % dxi, cosmo_grid % dyi, grid % clat,                     &
    1021           grid % clon, grid % ii, grid % jj, grid % w_horiz)
     998          grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz)
    1022999
    10231000!------------------------------------------------------------------------------
     
    10441021
    10451022       TYPE(grid_definition), POINTER    ::  grid
    1046        REAL ::  lonmin_pos,lonmax_pos, latmin_pos, latmax_pos
     1023       REAL(dp)                          ::  lonmin_pos,lonmax_pos, latmin_pos, latmax_pos
     1024       REAL(dp)                          ::  cosmo_dxi, cosmo_dyi
     1025
     1026       cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0))
     1027       cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0))
    10471028
    10481029       ! find horizontal index ranges for profile averaging
    1049        lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_grid % dxi
    1050        lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_grid % dxi
    1051        latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_grid % dyi
    1052        latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_grid % dyi
     1030       lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
     1031       lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
     1032       latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
     1033       latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
    10531034
    10541035       imin = FLOOR(lonmin_pos)
     
    10761057! Description:
    10771058! ------------
    1078 !> Helper function that computes horizontal domain extend in x or y direction
    1079 !> such that the centres of a boundary grid fall at -dx/2 or lx + dx/2.
    1080 !>
    1081 !> Input parameters:
    1082 !> -----------------
    1083 !> dxy : grid spacing in x or y direction
    1084 !> lxy : domain length in dxy direction
    1085 !>
    1086 !> Output parameters:
    1087 !> ------------------
    1088 !> boundary_extent : Domain minimum xymin (maximum xymax) if dxy < 0 (> 0)
    1089 !------------------------------------------------------------------------------!
    1090     REAL(dp) FUNCTION boundary_extent(dxy, lxy)
    1091         REAL(dp), INTENT(IN) ::  dxy, lxy
    1092 
    1093         boundary_extent = 0.5_dp * lxy + SIGN(lxy + ABS(dxy), dxy)
    1094 
    1095     END FUNCTION boundary_extent
    1096 
    1097 
    1098 !------------------------------------------------------------------------------!
    1099 ! Description:
    1100 ! ------------
    11011059!> Initializes grid_definition-type variables.
    11021060!>
    11031061!> Input parameters:
    11041062!> -----------------
    1105 !> mode : Initialization mode, distinguishes between PALM-4U and COSMO-DE grids
    1106 !>    as well as grids covering the boundary surfaces. Valid modes are:
     1063!> kind : Grid kind, distinguishes between PALM-4U and COSMO-DE grids
     1064!>    as well as grids covering the boundary surfaces. Valid kinds are:
    11071065!>       - 'palm'
    11081066!>       - 'cosmo-de'
     
    11141072!>    PALM-4U computational domain (i.e. the outer cell faces). The coordinates
    11151073!>    of the generated grid will be inferred from this information taking into
    1116 !>    account the initialization mode. For example, the coordinates of a
     1074!>    account the initialization mode ic_mode. For example, the coordinates of a
    11171075!>    boundary grid initialized using mode 'eastwest-scalar' will be located in
    11181076!>    planes one half grid point outwards of xmin and xmax.
     
    11261084!> grid : Grid variable to be initialized.
    11271085!------------------------------------------------------------------------------!
    1128     SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax, zmin, zmax,  &
    1129                                     x0, y0, z0, nx, ny, nz, dx, dy, dz, grid, mode)
     1086    SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax,              &
     1087                                    x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode)
    11301088        CHARACTER(LEN=*), INTENT(IN)           ::  kind
    1131         CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  mode
     1089        CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  ic_mode
    11321090        INTEGER, INTENT(IN)                    ::  nx, ny, nz
    1133         REAL(dp), INTENT(IN)                   ::  xmin, xmax, ymin, ymax, zmin, zmax
     1091        REAL(dp), INTENT(IN)                   ::  xmin, xmax, ymin, ymax
    11341092        REAL(dp), INTENT(IN)                   ::  x0, y0, z0
    1135         REAL(dp), OPTIONAL, INTENT(IN)         ::  dx, dy, dz
     1093        REAL(dp), INTENT(IN), TARGET, OPTIONAL ::  z(:)
     1094        REAL(dp), INTENT(IN), TARGET, OPTIONAL ::  zw(:)
    11361095        TYPE(grid_definition), INTENT(INOUT)   ::  grid
    11371096
     
    11421101        grid % lx = xmax - xmin
    11431102        grid % ly = ymax - ymin
    1144         grid % lz = zmax - zmin
    11451103
    11461104        grid % x0 = x0
     
    11511109
    11521110        CASE('boundary')
    1153            IF (.NOT. PRESENT(dx))  THEN
    1154               message = "dx is not present but needed for 'eastwest-scalar' "//&
    1155                         "initializaton."
     1111
     1112           IF (.NOT.PRESENT(z))  THEN
     1113              message = "z has not been passed but is required for 'boundary' grids"
    11561114              CALL abort('init_grid_definition', message)
    11571115           END IF
    1158            IF (.NOT. PRESENT(dy))  THEN
    1159               message = "dy is not present but needed for 'eastwest-scalar' "//&
    1160                         "initializaton."
    1161               CALL abort('init_grid_definition', message)
    1162            END IF
    1163            IF (.NOT. PRESENT(dz))  THEN
    1164               message = "dz is not present but needed for 'eastwest-scalar' "//&
    1165                         "initializaton."
    1166               CALL abort('init_grid_definition', message)
    1167            END IF
    1168 
    1169            grid % dx  = dx
    1170            grid % dy  = dy
    1171            grid % dz  = dz
    1172 
    1173            grid % dxi = 1.0_dp / grid % dx
    1174            grid % dyi = 1.0_dp / grid % dy
    1175            grid % dzi = 1.0_dp / grid % dz
    11761116
    11771117           ALLOCATE( grid % x(0:nx) )
     
    11811121           CALL linspace(ymin, ymax, grid % y)
    11821122
    1183            ALLOCATE( grid % z(0:nz) )
    1184            CALL linspace(zmin, zmax, grid % z)
     1123           grid % z => z
    11851124
    11861125           ! Allocate neighbour indices and weights
    1187            IF (TRIM(mode) .NE. 'profile')  THEN
    1188               ALLOCATE( grid % kk(0:nx, 0:ny, 0:nz, 2) )
     1126           IF (TRIM(ic_mode) .NE. 'profile')  THEN
     1127              ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
    11891128              grid % kk(:,:,:,:) = -1
    11901129
    1191               ALLOCATE( grid % w_verti(0:nx, 0:ny, 0:nz, 2) )
     1130              ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
    11921131              grid % w_verti(:,:,:,:) = 0.0_dp
    11931132           END IF
    11941133       
    11951134        CASE('boundary intermediate')
    1196            IF (.NOT. PRESENT(dx))  THEN
    1197               message = "dx is not present but needed for 'eastwest-scalar' "//&
    1198                         "initializaton."
    1199               CALL abort('init_grid_definition', message)
    1200            END IF
    1201            IF (.NOT. PRESENT(dy))  THEN
    1202               message = "dy is not present but needed for 'eastwest-scalar' "//&
    1203                         "initializaton."
    1204               CALL abort('init_grid_definition', message)
    1205            END IF
    1206            IF (.NOT. PRESENT(dz))  THEN
    1207               message = "dz is not present but needed for 'eastwest-scalar' "//&
    1208                         "initializaton."
    1209               CALL abort('init_grid_definition', message)
    1210            END IF
    1211 
    1212            grid % dx  = dx
    1213            grid % dy  = dy
    1214            grid % dz  = dz
    1215 
    1216            grid % dxi = 1.0_dp / grid % dx
    1217            grid % dyi = 1.0_dp / grid % dy
    1218            grid % dzi = 1.0_dp / grid % dz
    12191135
    12201136           ALLOCATE( grid % x(0:nx) )
     
    12241140           CALL linspace(ymin, ymax, grid % y)
    12251141
    1226            ALLOCATE( grid % z(0:nz) )
    1227            CALL linspace(zmin, zmax, grid % z)
    1228 
    12291142           ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny)  )
     1143
    12301144           CALL rotate_to_cosmo(                                               &
    1231               phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
    1232               lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude
     1145              phir = project( grid % y, y0, EARTH_RADIUS ) ,                   & ! = plate-carree latitude
     1146              lamr = project( grid % x, x0, EARTH_RADIUS ) ,                   & ! = plate-carree longitude
    12331147              phip = phi_cn, lamp = lambda_cn,                                 &
    12341148              phi  = grid % clat,                                              &
     
    12491163        ! corresponding latitudes and longitudes of the rotated pole grid.
    12501164        CASE('palm')
     1165
     1166           IF (.NOT.PRESENT(z))  THEN
     1167              message = "z has not been passed but is required for 'palm' grids"
     1168              CALL abort('init_grid_definition', message)
     1169           END IF
     1170
     1171           IF (.NOT.PRESENT(zw))  THEN
     1172              message = "zw has not been passed but is required for 'palm' grids"
     1173              CALL abort('init_grid_definition', message)
     1174           END IF
     1175
    12511176           grid % name(1) = 'x and lon'
    12521177           grid % name(2) = 'y and lat'
    12531178           grid % name(3) = 'z'
    12541179
    1255            grid % dx  = grid % lx / (nx + 1)
    1256            grid % dy  = grid % ly / (ny + 1)
    1257            grid % dz  = grid % lz / (nz + 1)
    1258 
    1259            grid % dxi = 1.0_dp / grid % dx
    1260            grid % dyi = 1.0_dp / grid % dy
    1261            grid % dzi = 1.0_dp / grid % dz
    1262 
    1263            ALLOCATE( grid % x(0:nx),   grid % y(0:ny),  grid % z(0:nz)  )
    1264            ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny), grid % zw(1:nz) )
    1265            CALL linspace(xmin + 0.5_dp*grid % dx, xmax - 0.5_dp*grid % dx, grid % x)
    1266            CALL linspace(ymin + 0.5_dp*grid % dy, ymax - 0.5_dp*grid % dy, grid % y)
    1267            CALL linspace(zmin + 0.5_dp*grid % dz, zmax - 0.5_dp*grid % dz, grid % z)
    1268            CALL linspace(xmin + grid % dx, xmax - grid % dx, grid % xu)
    1269            CALL linspace(ymin + grid % dy, ymax - grid % dy, grid % yv)
    1270            CALL linspace(zmin + grid % dz, zmax - grid % dz, grid % zw)
     1180           !TODO: Remove use of global dx, dy, dz variables. Consider
     1181           !TODO:   associating global x,y, and z arrays.
     1182           ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
     1183           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
     1184           CALL linspace(xmin + 0.5_dp* dx, xmax - 0.5_dp* dx, grid % x)
     1185           CALL linspace(ymin + 0.5_dp* dy, ymax - 0.5_dp* dy, grid % y)
     1186           grid % z => z
     1187           CALL linspace(xmin +  dx, xmax -  dx, grid % xu)
     1188           CALL linspace(ymin +  dy, ymax -  dy, grid % yv)
     1189           grid % zw => zw
    12711190
    12721191           grid % depths => depths
    12731192
    12741193           ! Allocate neighbour indices and weights
    1275            IF (TRIM(mode) .NE. 'profile')  THEN
    1276               ALLOCATE( grid % kk(0:nx, 0:ny, 0:nz, 2) )
     1194           IF (TRIM(ic_mode) .NE. 'profile')  THEN
     1195              ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
    12771196              grid % kk(:,:,:,:) = -1
    12781197
    1279               ALLOCATE( grid % w_verti(0:nx, 0:ny, 0:nz, 2) )
     1198              ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
    12801199              grid % w_verti(:,:,:,:) = 0.0_dp
    12811200           END IF
    12821201
    12831202        CASE('palm intermediate')
     1203
    12841204           grid % name(1) = 'x and lon'
    12851205           grid % name(2) = 'y and lat'
    1286            grid % name(3) = 'z'
    1287 
    1288            grid % dx  = grid % lx / (nx + 1)
    1289            grid % dy  = grid % ly / (ny + 1)
    1290            grid % dz  = grid % lz / (nz + 1)
    1291 
    1292            grid % dxi = 1.0_dp / grid % dx
    1293            grid % dyi = 1.0_dp / grid % dy
    1294            grid % dzi = 1.0_dp / grid % dz
    1295 
    1296            ALLOCATE( grid % x(0:nx),   grid % y(0:ny),  grid % z(0:nz)  )
    1297            ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny), grid % zw(1:nz) )
    1298            CALL linspace(xmin + 0.5_dp*grid % dx, xmax - 0.5_dp*grid % dx, grid % x)
    1299            CALL linspace(ymin + 0.5_dp*grid % dy, ymax - 0.5_dp*grid % dy, grid % y)
    1300            CALL linspace(zmin + 0.5_dp*grid % dz, zmax - 0.5_dp*grid % dz, grid % z)
    1301            CALL linspace(xmin + grid % dx, xmax - grid % dx, grid % xu)
    1302            CALL linspace(ymin + grid % dy, ymax - grid % dy, grid % yv)
    1303            CALL linspace(zmin + grid % dz, zmax - grid % dz, grid % zw)
     1206           grid % name(3) = 'interpolated hhl or hfl'
     1207
     1208           !TODO: Remove use of global dx, dy, dz variables. Consider
     1209           !TODO:   associating global x,y, and z arrays.
     1210           ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
     1211           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
     1212           CALL linspace(xmin + 0.5_dp*dx, xmax - 0.5_dp*dx, grid % x)
     1213           CALL linspace(ymin + 0.5_dp*dy, ymax - 0.5_dp*dy, grid % y)
     1214           CALL linspace(xmin + dx, xmax - dx, grid % xu)
     1215           CALL linspace(ymin + dy, ymax - dy, grid % yv)
    13041216
    13051217           grid % depths => depths
     
    13551267           grid % name(3) = 'height'
    13561268
    1357            grid % dx  = grid % lx / nx     ! = 0.025 deg, stored in radians
    1358            grid % dy  = grid % ly / ny     ! = 0.025 deg, stored in radians
    1359            grid % dz  = 0.0_dp             ! not defined yet
    1360 
    1361            grid % dxi = 1.0_dp / grid % dx ! [rad^-1]
    1362            grid % dyi = 1.0_dp / grid % dy ! [rad^-1]
    1363            grid % dzi = 0.0_dp             ! not defined yet
    1364 
    1365            ALLOCATE( grid % lon(0:nx),   grid % lat(0:ny),  grid % z(0:nz)  )
    1366            ALLOCATE( grid % lonu(0:nx),  grid % latv(0:ny), grid % zw(0:nz) )
     1269           ALLOCATE( grid % lon(0:nx),   grid % lat(0:ny)  )
     1270           ALLOCATE( grid % lonu(0:nx),  grid % latv(0:ny) )
    13671271
    13681272           CALL linspace(xmin, xmax, grid % lon)
    13691273           CALL linspace(ymin, ymax, grid % lat)
    1370            grid % lonu(:) = grid % lon + 0.5_dp * grid % dx
    1371            grid % latv(:) = grid % lat + 0.5_dp * grid % dy
     1274           grid % lonu(:) = grid % lon + 0.5_dp * (grid % lx / grid % nx)
     1275           grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny)
    13721276
    13731277           ! Point to heights of half levels (hhl) and compute heights of full
     
    13841288
    13851289    END SUBROUTINE init_grid_definition
     1290
     1291
     1292!------------------------------------------------------------------------------!
     1293! Description:
     1294! ------------
     1295!> PALM's stretched vertical grid generator. Forked from PALM revision 3139, see
     1296!> https://palm.muk.uni-hannover.de/trac/browser/palm/trunk/SOURCE/init_grid.f90?rev=3139
     1297!>
     1298!> This routine computes the levels of scalar points. The levels of the velocity
     1299!> points are then obtained as the midpoints inbetween using the INIFOR routine
     1300!> 'modpoints'.
     1301!------------------------------------------------------------------------------!
     1302    SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, &
     1303                           dz_stretch_level_start, dz_stretch_level_end,       &
     1304                           dz_stretch_factor_array)
     1305
     1306       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  z, dz, dz_stretch_factor_array
     1307       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_level_start, dz_stretch_level_end
     1308       REAL(dp), INTENT(IN) ::  dz_max, dz_stretch_factor, dz_stretch_level
     1309
     1310       INTEGER ::  number_stretch_level_start        !< number of user-specified start levels for stretching
     1311       INTEGER ::  number_stretch_level_end          !< number of user-specified end levels for stretching
     1312
     1313       REAL(dp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end
     1314       REAL(dp) ::  dz_level_end, dz_stretched
     1315
     1316       INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
     1317       INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
     1318       INTEGER ::  dz_stretch_level_index = 0
     1319       INTEGER ::  k, n, number_dz
     1320!
     1321!-- Compute height of u-levels from constant grid length and dz stretch factors
     1322       IF ( dz(1) == -1.0_dp )  THEN
     1323          message = 'missing dz'
     1324          CALL abort( 'stretched_z', message)
     1325       ELSEIF ( dz(1) <= 0.0_dp )  THEN
     1326          WRITE( message, * ) 'dz=', dz(1),' <= 0.0'
     1327          CALL abort( 'stretched_z', message)
     1328       ENDIF
     1329
     1330!
     1331!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
     1332!-- if it was set by the user
     1333       IF ( dz_stretch_level /= -9999999.9_dp ) THEN
     1334          dz_stretch_level_start(1) = dz_stretch_level
     1335       ENDIF
     1336       
     1337!
     1338!-- Determine number of dz values and stretching levels specified by the
     1339!-- user to allow right controlling of the stretching mechanism and to
     1340!-- perform error checks. The additional requirement that dz /= dz_max
     1341!-- for counting number of user-specified dz values is necessary. Otherwise
     1342!-- restarts would abort if the old stretching mechanism with dz_stretch_level
     1343!-- is used (Attention: The user is not allowed to specify a dz value equal
     1344!-- to the default of dz_max = 999.0).
     1345       number_dz = COUNT( dz /= -1.0_dp .AND. dz /= dz_max )
     1346       number_stretch_level_start = COUNT( dz_stretch_level_start /=           &
     1347                                           -9999999.9_dp )
     1348       number_stretch_level_end = COUNT( dz_stretch_level_end /=               &
     1349                                         9999999.9_dp )
     1350
     1351!
     1352!-- The number of specified end levels +1 has to be the same than the number
     1353!-- of specified dz values
     1354       IF ( number_dz /= number_stretch_level_end + 1 ) THEN
     1355          WRITE( message, * ) 'The number of values for dz = ',                &
     1356                              number_dz, 'has to be the same than ',           &
     1357                              'the number of values for ',                     &
     1358                              'dz_stretch_level_end + 1 = ',                   &
     1359                              number_stretch_level_end+1
     1360          CALL abort( 'stretched_z', message)
     1361       ENDIF
     1362   
     1363!
     1364!--    The number of specified start levels has to be the same or one less than
     1365!--    the number of specified dz values
     1366       IF ( number_dz /= number_stretch_level_start + 1 .AND.                  &
     1367            number_dz /= number_stretch_level_start ) THEN
     1368          WRITE( message, * ) 'The number of values for dz = ',         &
     1369                              number_dz, 'has to be the same or one ', &
     1370                              'more than& the number of values for ',  &
     1371                              'dz_stretch_level_start = ',             &
     1372                              number_stretch_level_start
     1373          CALL abort( 'stretched_z', message)
     1374       ENDIF
     1375   
     1376!--    The number of specified start levels has to be the same or one more than
     1377!--    the number of specified end levels
     1378       IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.   &
     1379            number_stretch_level_start /= number_stretch_level_end ) THEN
     1380          WRITE( message, * ) 'The number of values for ',              &
     1381                              'dz_stretch_level_start = ',              &
     1382                              dz_stretch_level_start, 'has to be the ',&
     1383                              'same or one more than& the number of ', &
     1384                              'values for dz_stretch_level_end = ',    &
     1385                              number_stretch_level_end
     1386          CALL abort( 'stretched_z', message)
     1387       ENDIF
     1388
     1389!
     1390!-- Initialize dz for the free atmosphere with the value of dz_max
     1391       IF ( dz(number_stretch_level_start+1) == -1.0_dp .AND.                     &
     1392            number_stretch_level_start /= 0 ) THEN
     1393          dz(number_stretch_level_start+1) = dz_max
     1394       ENDIF
     1395       
     1396!
     1397!-- Initialize the stretching factor if (infinitely) stretching in the free
     1398!-- atmosphere is desired (dz_stretch_level_end was not specified for the
     1399!-- free atmosphere)
     1400       IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN
     1401          dz_stretch_factor_array(number_stretch_level_start) =                   &
     1402          dz_stretch_factor
     1403       ENDIF
     1404
     1405!-- Allocation of arrays for stretching
     1406       ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
     1407
     1408!
     1409!--    The stretching region has to be large enough to allow for a smooth
     1410!--    transition between two different grid spacings
     1411       DO n = 1, number_stretch_level_start
     1412          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
     1413                                        4 * MAX( dz(n),dz(n+1) )
     1414       ENDDO
     1415
     1416       IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
     1417                 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
     1418       !IF ( ANY( min_dz_stretch_level_end >      &
     1419       !          dz_stretch_level_end ) ) THEN
     1420             message = 'Each dz_stretch_level_end has to be larger '  // &
     1421                       'than its corresponding value for ' //            &
     1422                       'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
     1423                       'to allow for smooth grid stretching'
     1424             CALL abort('stretched_z', message)
     1425       ENDIF
     1426       
     1427!
     1428!--    Stretching must not be applied within the prandtl_layer
     1429!--    (first two grid points). For the default case dz_stretch_level_start
     1430!--    is negative. Therefore the absolut value is checked here.
     1431       IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_dp ) ) THEN
     1432          WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',&
     1433                              'larger than ', dz(1) * 1.5
     1434             CALL abort( 'stretched_z', message)
     1435       ENDIF
     1436
     1437!
     1438!--    The stretching has to start and end on a grid level. Therefore
     1439!--    user-specified values have to ''interpolate'' to the next lowest level
     1440       IF ( number_stretch_level_start /= 0 ) THEN
     1441          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
     1442                                            dz(1)/2.0) / dz(1) )               &
     1443                                      * dz(1) + dz(1)/2.0
     1444       ENDIF
     1445       
     1446       IF ( number_stretch_level_start > 1 ) THEN
     1447          DO n = 2, number_stretch_level_start
     1448             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
     1449                                              dz(n) ) * dz(n)
     1450          ENDDO
     1451       ENDIF
     1452       
     1453       IF ( number_stretch_level_end /= 0 ) THEN
     1454          DO n = 1, number_stretch_level_end
     1455             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
     1456                                            dz(n+1) ) * dz(n+1)
     1457          ENDDO
     1458       ENDIF
     1459 
     1460!
     1461!--    Determine stretching factor if necessary
     1462       IF ( number_stretch_level_end >= 1 ) THEN
     1463          CALL calculate_stretching_factor( number_stretch_level_end, dz,      &
     1464                                            dz_stretch_factor,                 &
     1465                                            dz_stretch_factor_array,           &   
     1466                                            dz_stretch_level_end,              &
     1467                                            dz_stretch_level_start )
     1468       ENDIF
     1469
     1470       z(1) = dz(1) * 0.5_dp
     1471!
     1472       dz_stretch_level_index = n
     1473       dz_stretched = dz(1)
     1474       DO  k = 2, n
     1475
     1476          IF ( dz_stretch_level <= z(k-1)  .AND.  dz_stretched < dz_max )  THEN
     1477
     1478             dz_stretched = dz_stretched * dz_stretch_factor
     1479             dz_stretched = MIN( dz_stretched, dz_max )
     1480
     1481             IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1
     1482
     1483          ENDIF
     1484
     1485          z(k) = z(k-1) + dz_stretched
     1486
     1487       ENDDO
     1488!--    Determine u and v height levels considering the possibility of grid
     1489!--    stretching in several heights.
     1490       n = 1
     1491       dz_stretch_level_start_index(:) = UBOUND(z, 1)
     1492       dz_stretch_level_end_index(:) = UBOUND(z, 1)
     1493       dz_stretched = dz(1)
     1494
     1495!--    The default value of dz_stretch_level_start is negative, thus the first
     1496!--    condition is always true. Hence, the second condition is necessary.
     1497       DO  k = 2, UBOUND(z, 1)
     1498          IF ( dz_stretch_level_start(n) <= z(k-1) .AND.                      &
     1499               dz_stretch_level_start(n) /= -9999999.9_dp ) THEN
     1500             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
     1501             
     1502             IF ( dz(n) > dz(n+1) ) THEN
     1503                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
     1504             ELSE
     1505                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
     1506             ENDIF
     1507             
     1508             IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) )            &
     1509             dz_stretch_level_start_index(n) = k-1
     1510             
     1511          ENDIF
     1512         
     1513          z(k) = z(k-1) + dz_stretched
     1514         
     1515!
     1516!--       Make sure that the stretching ends exactly at dz_stretch_level_end
     1517          dz_level_end = ABS( z(k) - dz_stretch_level_end(n) )
     1518         
     1519          IF ( dz_level_end < dz(n+1)/3.0 ) THEN
     1520             z(k) = dz_stretch_level_end(n)
     1521             dz_stretched = dz(n+1)
     1522             dz_stretch_level_end_index(n) = k
     1523             n = n + 1             
     1524          ENDIF
     1525       ENDDO
     1526
     1527       DEALLOCATE( min_dz_stretch_level_end )
     1528
     1529    END SUBROUTINE stretched_z
     1530
     1531
     1532! Description: [PALM subroutine]
     1533! -----------------------------------------------------------------------------!
     1534!> Calculation of the stretching factor through an iterative method. Ideas were
     1535!> taken from the paper "Regional stretched grid generation and its application
     1536!> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
     1537!> system of equations has two variables (r,l) but four requirements
     1538!> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
     1539!> results into an overdetermined system.
     1540!------------------------------------------------------------------------------!
     1541 SUBROUTINE calculate_stretching_factor( number_end, dz, dz_stretch_factor,    &
     1542                                         dz_stretch_factor_array,              &   
     1543                                         dz_stretch_level_end,                 &
     1544                                         dz_stretch_level_start )
     1545 
     1546    REAL(dp), DIMENSION(:), INTENT(IN)    ::  dz
     1547    REAL(dp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_factor_array
     1548    REAL(dp), DIMENSION(:), INTENT(IN)    ::  dz_stretch_level_end, dz_stretch_level_start
     1549    REAL(dp)                              ::  dz_stretch_factor
     1550 
     1551    INTEGER ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached 
     1552    INTEGER ::  l_rounded   !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2
     1553    INTEGER ::  n           !< loop variable for stretching
     1554   
     1555    INTEGER, INTENT(IN) ::  number_end !< number of user-specified end levels for stretching
     1556       
     1557    REAL(dp) ::  delta_l               !< absolute difference between l and l_rounded
     1558    REAL(dp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
     1559    REAL(dp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
     1560    REAL(dp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration
     1561    REAL(dp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
     1562    REAL(dp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
     1563    REAL(dp) ::  numerator             !< numerator of the quotient
     1564    REAL(dp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
     1565    REAL(dp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
     1566   
     1567    REAL(dp) ::  dz_stretch_factor_array_2(9) = 1.08_dp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
     1568   
     1569    REAL(dp), PARAMETER ::  stretch_factor_interval = 1.0E-06  !< interval for sampling possible stretching factors
     1570    REAL(dp), PARAMETER ::  stretch_factor_lower_limit = 0.88  !< lowest possible stretching factor
     1571    REAL(dp), PARAMETER ::  stretch_factor_upper_limit = 1.12  !< highest possible stretching factor
     1572 
     1573 
     1574    l = 0
     1575    DO  n = 1, number_end
     1576   
     1577       iterations = 1
     1578       stretch_factor_1 = 1.0
     1579       stretch_factor_2 = 1.0
     1580       delta_total_old = 1.0
     1581       
     1582       IF ( dz(n) > dz(n+1) ) THEN
     1583          DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )
     1584             
     1585             stretch_factor_1 = 1.0 - iterations * stretch_factor_interval
     1586             distance = ABS( dz_stretch_level_end(n) -                   &
     1587                        dz_stretch_level_start(n) )
     1588             numerator = distance*stretch_factor_1/dz(n) +               &
     1589                         stretch_factor_1 - distance/dz(n)
     1590             
     1591             IF ( numerator > 0.0 ) THEN
     1592                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
     1593                l_rounded = NINT( l )
     1594                delta_l = ABS( l_rounded - l ) / l
     1595             ENDIF
     1596             
     1597             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
     1598             
     1599             delta_stretch_factor = ABS( stretch_factor_1 -              &
     1600                                         stretch_factor_2 ) /            &
     1601                                    stretch_factor_2
     1602             
     1603             delta_total_new = delta_l + delta_stretch_factor
     1604
     1605!
     1606!--                stretch_factor_1 is taken to guarantee that the stretching
     1607!--                procedure ends as close as possible to dz_stretch_level_end.
     1608!--                stretch_factor_2 would guarantee that the stretched dz(n) is
     1609!--                equal to dz(n+1) after l_rounded grid levels.
     1610             IF (delta_total_new < delta_total_old) THEN
     1611                dz_stretch_factor_array(n) = stretch_factor_1
     1612                dz_stretch_factor_array_2(n) = stretch_factor_2
     1613                delta_total_old = delta_total_new
     1614             ENDIF
     1615             
     1616             iterations = iterations + 1
     1617           
     1618          ENDDO
     1619             
     1620       ELSEIF ( dz(n) < dz(n+1) ) THEN
     1621          DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
     1622                     
     1623             stretch_factor_1 = 1.0 + iterations * stretch_factor_interval
     1624             distance = ABS( dz_stretch_level_end(n) -                      &
     1625                        dz_stretch_level_start(n) )
     1626             numerator = distance*stretch_factor_1/dz(n) +                  &
     1627                         stretch_factor_1 - distance/dz(n)
     1628             
     1629             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
     1630             l_rounded = NINT( l )
     1631             delta_l = ABS( l_rounded - l ) / l
     1632             
     1633             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
     1634
     1635             delta_stretch_factor = ABS( stretch_factor_1 -                 &
     1636                                        stretch_factor_2 ) /                &
     1637                                        stretch_factor_2
     1638             
     1639             delta_total_new = delta_l + delta_stretch_factor
     1640             
     1641!
     1642!--                stretch_factor_1 is taken to guarantee that the stretching
     1643!--                procedure ends as close as possible to dz_stretch_level_end.
     1644!--                stretch_factor_2 would guarantee that the stretched dz(n) is
     1645!--                equal to dz(n+1) after l_rounded grid levels.
     1646             IF (delta_total_new < delta_total_old) THEN
     1647                dz_stretch_factor_array(n) = stretch_factor_1
     1648                dz_stretch_factor_array_2(n) = stretch_factor_2
     1649                delta_total_old = delta_total_new
     1650             ENDIF
     1651             
     1652             iterations = iterations + 1
     1653          ENDDO
     1654         
     1655       ELSE
     1656          message = 'Two adjacent values of dz must be different'
     1657          CALL abort( 'calculate_stretching_factor', message)
     1658       ENDIF
     1659
     1660!
     1661!--    Check if also the second stretching factor fits into the allowed
     1662!--    interval. If not, print a warning for the user.
     1663       IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR.     &
     1664            dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN
     1665          WRITE( message, * ) 'stretch_factor_2 = ',                    &
     1666                                     dz_stretch_factor_array_2(n), ' which is',&
     1667                                     ' responsible for exactly reaching& dz =',&
     1668                                      dz(n+1), 'after a specific amount of',   &
     1669                                     ' grid levels& exceeds the upper',        &
     1670                                     ' limit =', stretch_factor_upper_limit,   &
     1671                                     ' &or lower limit = ',                    &
     1672                                     stretch_factor_lower_limit
     1673          CALL abort( 'calculate_stretching_factor', message )
     1674           
     1675       ENDIF
     1676    ENDDO
     1677       
     1678 END SUBROUTINE calculate_stretching_factor
     1679
     1680    SUBROUTINE midpoints(z, zw)
     1681
     1682        REAL(dp), INTENT(IN)  ::  z(0:)
     1683        REAL(dp), INTENT(OUT) ::  zw(1:)
     1684
     1685        INTEGER ::  k
     1686
     1687        DO k = 1, UBOUND(zw, 1)
     1688           zw(k) = 0.5_dp * (z(k-1) + z(k))
     1689        END DO
     1690
     1691    END SUBROUTINE midpoints
    13861692
    13871693
     
    14091715       )
    14101716
    1411        !potential temperature, surface pressure
     1717       !potential temperature, surface pressure, including nudging and subsidence
    14121718       io_group_list(3) = init_io_group(                                       &
    14131719          in_files = flow_files,                                               &
    1414           out_vars = [output_var_table(3:8), output_var_table(42:42)],         &
     1720          out_vars = [output_var_table(3:8), output_var_table(42:42),          &
     1721                      output_var_table(49:51)],                                &
    14151722          in_var_list = (/input_var_table(3), input_var_table(17)/),           &
    14161723          kind = 'temperature'                                                 &
    14171724       )
    14181725
    1419        !specific humidity
     1726       !specific humidity including nudging and subsidence
    14201727       io_group_list(4) = init_io_group(                                       &
    14211728          in_files = flow_files,                                               &
    1422           out_vars = output_var_table(9:14),                                   &
     1729          out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
    14231730          in_var_list = input_var_table(4:4),                                  &
    14241731          kind = 'scalar'                                                      &
     
    14281735       io_group_list(5) = init_io_group(                                       &
    14291736          in_files = flow_files,                                               &
    1430           out_vars = [output_var_table(15:26), output_var_table(43:44)],       &
     1737          out_vars = [output_var_table(15:26), output_var_table(43:46)],       &
    14311738          !out_vars = output_var_table(15:20),                                  &
    14321739          in_var_list = input_var_table(5:6),                                  &
     
    14441751       !io_group_list(6) % to_be_processed = .FALSE.
    14451752   
    1446        !w velocity
     1753       !w velocity and subsidence and w nudging
    14471754       io_group_list(7) = init_io_group(                                       &
    14481755          in_files = flow_files,                                               &
    1449           out_vars = output_var_table(27:32),                                  &
     1756          out_vars = [output_var_table(27:32), output_var_table(47:48)],       &
    14501757          in_var_list = input_var_table(7:7),                                  &
    14511758          kind = 'scalar'                                                      &
     
    14591766          kind = 'accumulated'                                                 &
    14601767       )
     1768       io_group_list(8) % to_be_processed = .FALSE.
    14611769
    14621770       !snow
     
    14671775          kind = 'accumulated'                                                 &
    14681776       )
     1777       io_group_list(9) % to_be_processed = .FALSE.
    14691778
    14701779       !graupel
     
    14751784          kind = 'accumulated'                                                 &
    14761785       )
     1786       io_group_list(10) % to_be_processed = .FALSE.
    14771787
    14781788       !evapotranspiration
     
    14831793          kind = 'accumulated'                                                 &
    14841794       )
     1795       io_group_list(11) % to_be_processed = .FALSE.
    14851796
    14861797       !2m air temperature
     
    15171828          kind = 'running average'                                             &
    15181829       )
     1830       io_group_list(15) % to_be_processed = .FALSE.
    15191831
    15201832       !lw radiation balance
     
    15251837          kind = 'running average'                                             &
    15261838       )
     1839       io_group_list(16) % to_be_processed = .FALSE.
    15271840
    15281841    END SUBROUTINE setup_io_groups
     
    15621875       CALL report('fini_grids', 'Deallocating grids')
    15631876       
     1877       DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
     1878
    15641879       DEALLOCATE(palm_grid%x,  palm_grid%y,  palm_grid%z,                     &
    15651880                  palm_grid%xu, palm_grid%yv, palm_grid%zw,                    &
     
    15721887                  palm_intermediate%clonu, palm_intermediate%clatu)
    15731888
    1574        DEALLOCATE(cosmo_grid%lon,  cosmo_grid%lat,  cosmo_grid%z,              &
    1575                   cosmo_grid%lonu, cosmo_grid%latv, cosmo_grid%zw,             &
     1889       DEALLOCATE(cosmo_grid%lon,  cosmo_grid%lat,                             &
     1890                  cosmo_grid%lonu, cosmo_grid%latv,                            &
    15761891                  cosmo_grid%hfl)
    15771892
     
    15841899!> Initializes the the variable list.
    15851900!------------------------------------------------------------------------------!
    1586     SUBROUTINE setup_variable_tables(mode)
    1587        CHARACTER(LEN=*), INTENT(IN) ::  mode
     1901    SUBROUTINE setup_variable_tables(ic_mode)
     1902       CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
    15881903       TYPE(nc_var), POINTER        ::  var
    15891904
    1590        IF (TRIM(start_date) == '')  THEN
     1905       IF (TRIM(cfg % start_date) == '')  THEN
    15911906          message = 'Simulation start date has not been set.'
    15921907          CALL abort('setup_variable_tables', message)
    15931908       END IF
    15941909
    1595        nc_source_text = 'COSMO-DE analysis from ' // TRIM(start_date)
     1910       nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg % start_date)
    15961911
    15971912       n_invar = 17
    1598        n_outvar = 44
     1913       n_outvar = 55
    15991914       ALLOCATE( input_var_table(n_invar) )
    16001915       ALLOCATE( output_var_table(n_outvar) )
     
    16932008!- Section 2: NetCDF output variables
    16942009!------------------------------------------------------------------------------
     2010!
     2011!------------------------------------------------------------------------------
     2012! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions
     2013!------------------------------------------------------------------------------
    16952014       output_var_table(1) = init_nc_var(                                      &
    16962015          name              = 'init_soil_t',                                   &
     
    17182037
    17192038       output_var_table(3) = init_nc_var(                                      &
    1720           name              = 'init_pt',                                       &
     2039          name              = 'init_atmosphere_pt',                            &
    17212040          std_name          = "",                                              &
    17222041          long_name         = "initial potential temperature",                 &
     
    17272046          grid              = palm_grid,                                       &
    17282047          intermediate_grid = palm_intermediate,                               &
    1729           is_profile = (TRIM(mode) == 'profile')                               &
    1730        )
    1731        IF (TRIM(mode) == 'profile')  THEN
     2048          is_profile = (TRIM(ic_mode) == 'profile')                            &
     2049       )
     2050       IF (TRIM(ic_mode) == 'profile')  THEN
    17322051          output_var_table(3) % grid => scalar_profile_grid
    17332052          output_var_table(3) % intermediate_grid => scalar_profile_intermediate
     
    17952114
    17962115       output_var_table(9) = init_nc_var(                                      &
    1797           name              = 'init_qv',                                       &
     2116          name              = 'init_atmosphere_qv',                            &
    17982117          std_name          = "",                                              &
    17992118          long_name         = "initial specific humidity",                     &
     
    18042123          grid              = palm_grid,                                       &
    18052124          intermediate_grid = palm_intermediate,                               &
    1806           is_profile = (TRIM(mode) == 'profile')                               &
    1807        )
    1808        IF (TRIM(mode) == 'profile')  THEN
     2125          is_profile = (TRIM(ic_mode) == 'profile')                            &
     2126       )
     2127       IF (TRIM(ic_mode) == 'profile')  THEN
    18092128          output_var_table(9) % grid => scalar_profile_grid
    18102129          output_var_table(9) % intermediate_grid => scalar_profile_intermediate
     
    18722191
    18732192       output_var_table(15) = init_nc_var(                                     &
    1874           name              = 'init_u',                                        &
     2193          name              = 'init_atmosphere_u',                             &
    18752194          std_name          = "",                                              &
    18762195          long_name         = "initial wind component in x direction",         &
     
    18812200          grid              = u_initial_grid,                                  &
    18822201          intermediate_grid = u_initial_intermediate,                          &
    1883           is_profile = (TRIM(mode) == 'profile')                               &
    1884        )
    1885        IF (TRIM(mode) == 'profile')  THEN
     2202          is_profile = (TRIM(ic_mode) == 'profile')                            &
     2203       )
     2204       IF (TRIM(ic_mode) == 'profile')  THEN
    18862205          output_var_table(15) % grid => scalar_profile_grid
    18872206          output_var_table(15) % intermediate_grid => scalar_profile_intermediate
     
    19492268
    19502269       output_var_table(21) = init_nc_var(                                     &
    1951           name              = 'init_v',                                        &
     2270          name              = 'init_atmosphere_v',                             &
    19522271          std_name          = "",                                              &
    19532272          long_name         = "initial wind component in y direction",         &
     
    19582277          grid              = v_initial_grid,                                  &
    19592278          intermediate_grid = v_initial_intermediate,                          &
    1960           is_profile = (TRIM(mode) == 'profile')                               &
    1961        )
    1962        IF (TRIM(mode) == 'profile')  THEN
     2279          is_profile = (TRIM(ic_mode) == 'profile')                            &
     2280       )
     2281       IF (TRIM(ic_mode) == 'profile')  THEN
    19632282          output_var_table(21) % grid => scalar_profile_grid
    19642283          output_var_table(21) % intermediate_grid => scalar_profile_intermediate
     
    20262345
    20272346       output_var_table(27) = init_nc_var(                                     &
    2028           name              = 'init_w',                                        &
     2347          name              = 'init_atmosphere_w',                             &
    20292348          std_name          = "",                                              &
    20302349          long_name         = "initial wind component in z direction",         &
     
    20352354          grid              = w_initial_grid,                                  &
    20362355          intermediate_grid = w_initial_intermediate,                          &
    2037           is_profile = (TRIM(mode) == 'profile')                               &
    2038        )
    2039        IF (TRIM(mode) == 'profile')  THEN
     2356          is_profile = (TRIM(ic_mode) == 'profile')                            &
     2357       )
     2358       IF (TRIM(ic_mode) == 'profile')  THEN
    20402359          output_var_table(27) % grid => w_profile_grid
    20412360          output_var_table(27) % intermediate_grid => w_profile_intermediate
     
    22092528          intermediate_grid = palm_intermediate                                &
    22102529       )
    2211 
     2530!
     2531!------------------------------------------------------------------------------
     2532! Section 2.2: Idealized large-scale forcings
     2533!------------------------------------------------------------------------------
    22122534       output_var_table(42) = init_nc_var(                                     &
    22132535          name              = 'surface_forcing_surface_pressure',              &
     
    22272549          long_name         = "geostrophic wind (u component)",                &
    22282550          units             = "m/s",                                           &
    2229           kind              = "profile",                                       &
    2230           input_id          = 1,                                               &
    2231           output_file       = output_file,                                     &
    2232           grid              = palm_grid,                                       &
    2233           intermediate_grid = palm_intermediate                                &
     2551          kind              = "constant scalar profile",                       &
     2552          input_id          = 1,                                               &
     2553          output_file       = output_file,                                     &
     2554          grid              = scalar_profile_grid,                             &
     2555          intermediate_grid = scalar_profile_intermediate                      &
    22342556       )
    22352557
     
    22392561          long_name         = "geostrophic wind (v component)",                &
    22402562          units             = "m/s",                                           &
    2241           kind              = "profile",                                       &
    2242           input_id          = 1,                                               &
    2243           output_file       = output_file,                                     &
    2244           grid              = palm_grid,                                       &
    2245           intermediate_grid = palm_intermediate                                &
    2246        )
     2563          kind              = "constant scalar profile",                       &
     2564          input_id          = 1,                                               &
     2565          output_file       = output_file,                                     &
     2566          grid              = scalar_profile_grid,                             &
     2567          intermediate_grid = scalar_profile_intermediate                      &
     2568       )
     2569
     2570       output_var_table(45) = init_nc_var(                                     &
     2571          name              = 'nudging_u',                                     &
     2572          std_name          = "",                                              &
     2573          long_name         = "wind component in x direction",                 &
     2574          units             = "m/s",                                           &
     2575          kind              = "large-scale scalar forcing",                    &
     2576          input_id          = 1,                                               &
     2577          output_file       = output_file,                                     &
     2578          grid              = scalar_profile_grid,                             &
     2579          intermediate_grid = scalar_profile_intermediate                      &
     2580       )
     2581       output_var_table(45) % to_be_processed = ls_forcing_variables_required
     2582
     2583       output_var_table(46) = init_nc_var(                                     &
     2584          name              = 'nudging_v',                                     &
     2585          std_name          = "",                                              &
     2586          long_name         = "wind component in y direction",                 &
     2587          units             = "m/s",                                           &
     2588          kind              = "large-scale scalar forcing",                    &
     2589          input_id          = 1,                                               &
     2590          output_file       = output_file,                                     &
     2591          grid              = scalar_profile_grid,                             &
     2592          intermediate_grid = scalar_profile_intermediate                      &
     2593       )
     2594       output_var_table(46) % to_be_processed = ls_forcing_variables_required
     2595
     2596       output_var_table(47) = init_nc_var(                                     &
     2597          name              = 'ls_forcing_sub_w',                              &
     2598          std_name          = "",                                              &
     2599          long_name         = "subsidence velocity of w",                      &
     2600          units             = "m/s",                                           &
     2601          kind              = "large-scale w forcing",                         &
     2602          input_id          = 1,                                               &
     2603          output_file       = output_file,                                     &
     2604          grid              = w_profile_grid,                                  &
     2605          intermediate_grid = w_profile_intermediate                           &
     2606       )
     2607       output_var_table(47) % to_be_processed = ls_forcing_variables_required
     2608
     2609       output_var_table(48) = init_nc_var(                                     &
     2610          name              = 'nudging_w',                                     &
     2611          std_name          = "",                                              &
     2612          long_name         = "wind component in w direction",                 &
     2613          units             = "m/s",                                           &
     2614          kind              = "large-scale w forcing",                         &
     2615          input_id          = 1,                                               &
     2616          output_file       = output_file,                                     &
     2617          grid              = w_profile_grid,                                  &
     2618          intermediate_grid = w_profile_intermediate                           &
     2619       )
     2620       output_var_table(48) % to_be_processed = ls_forcing_variables_required
     2621
     2622
     2623       output_var_table(49) = init_nc_var(                                     &
     2624          name              = 'ls_forcing_adv_pt',                             &
     2625          std_name          = "",                                              &
     2626          long_name         = "advection of potential temperature",            &
     2627          units             = "K/s",                                           &
     2628          kind              = "large-scale scalar forcing",                    &
     2629          input_id          = 1,                                               &
     2630          output_file       = output_file,                                     &
     2631          grid              = scalar_profile_grid,                             &
     2632          intermediate_grid = scalar_profile_intermediate                      &
     2633       )
     2634       output_var_table(49) % to_be_processed = ls_forcing_variables_required
     2635
     2636       output_var_table(50) = init_nc_var(                                     &
     2637          name              = 'ls_forcing_sub_pt',                             &
     2638          std_name          = "",                                              &
     2639          long_name         = "subsidence velocity of potential temperature",  &
     2640          units             = "K/s",                                           &
     2641          kind              = "large-scale scalar forcing",                    &
     2642          input_id          = 1,                                               &
     2643          output_file       = output_file,                                     &
     2644          grid              = scalar_profile_grid,                             &
     2645          intermediate_grid = scalar_profile_intermediate                      &
     2646       )
     2647       output_var_table(50) % to_be_processed = ls_forcing_variables_required
     2648
     2649       output_var_table(51) = init_nc_var(                                     &
     2650          name              = 'nudging_pt',                                    &
     2651          std_name          = "",                                              &
     2652          long_name         = "potential temperature",                         &
     2653          units             = "K",                                             &
     2654          kind              = "large-scale scalar forcing",                    &
     2655          input_id          = 1,                                               &
     2656          output_file       = output_file,                                     &
     2657          grid              = scalar_profile_grid,                             &
     2658          intermediate_grid = scalar_profile_intermediate                      &
     2659       )
     2660       output_var_table(51) % to_be_processed = ls_forcing_variables_required
     2661
     2662       output_var_table(52) = init_nc_var(                                     &
     2663          name              = 'ls_forcing_adv_qv',                             &
     2664          std_name          = "",                                              &
     2665          long_name         = "advection of specific humidity",                &
     2666          units             = "kg/kg/s",                                       &
     2667          kind              = "large-scale scalar forcing",                    &
     2668          input_id          = 1,                                               &
     2669          output_file       = output_file,                                     &
     2670          grid              = scalar_profile_grid,                             &
     2671          intermediate_grid = scalar_profile_intermediate                      &
     2672       )
     2673       output_var_table(52) % to_be_processed = ls_forcing_variables_required
     2674
     2675
     2676       output_var_table(53) = init_nc_var(                                     &
     2677          name              = 'ls_forcing_sub_qv',                             &
     2678          std_name          = "",                                              &
     2679          long_name         = "subsidence velocity of specific humidity",      &
     2680          units             = "kg/kg/s",                                       &
     2681          kind              = "large-scale scalar forcing",                    &
     2682          input_id          = 1,                                               &
     2683          output_file       = output_file,                                     &
     2684          grid              = scalar_profile_grid,                             &
     2685          intermediate_grid = scalar_profile_intermediate                      &
     2686       )
     2687       output_var_table(53) % to_be_processed = ls_forcing_variables_required
     2688
     2689       output_var_table(54) = init_nc_var(                                     &
     2690          name              = 'nudging_qv',                                    &
     2691          std_name          = "",                                              &
     2692          long_name         = "specific humidity",                             &
     2693          units             = "kg/kg",                                         &
     2694          kind              = "large-scale scalar forcing",                    &
     2695          input_id          = 1,                                               &
     2696          output_file       = output_file,                                     &
     2697          grid              = scalar_profile_grid,                             &
     2698          intermediate_grid = scalar_profile_intermediate                      &
     2699       )
     2700       output_var_table(54) % to_be_processed = ls_forcing_variables_required
     2701
     2702       output_var_table(55) = init_nc_var(                                     &
     2703          name              = 'nudging_tau',                                   &
     2704          std_name          = "",                                              &
     2705          long_name         = "nudging relaxation time scale",                 &
     2706          units             = "s",                                             &
     2707          kind              = "constant scalar profile",                       &
     2708          input_id          = 1,                                               &
     2709          output_file       = output_file,                                     &
     2710          grid              = scalar_profile_grid,                             &
     2711          intermediate_grid = scalar_profile_intermediate                      &
     2712       )
     2713       output_var_table(55) % to_be_processed = ls_forcing_variables_required
     2714
    22472715
    22482716       ! Attributes shared among all variables
    22492717       output_var_table(:) % source = nc_source_text
     2718
    22502719
    22512720    END SUBROUTINE setup_variable_tables
     
    22602729!------------------------------------------------------------------------------!
    22612730    FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id,     &
    2262                          grid, intermediate_grid, output_file, is_profile) RESULT(var)
     2731                         grid, intermediate_grid, output_file, is_profile      &
     2732             ) RESULT(var)
    22632733
    22642734       CHARACTER(LEN=*), INTENT(IN)      ::  name, std_name, long_name, units, kind
     
    23672837          var % task            = "average profile"
    23682838
    2369        CASE ( 'surface forcing' )
    2370           var % lod             = 1
     2839       CASE( 'surface forcing' )
     2840          var % lod             = -1
    23712841          var % ndim            = 3
    23722842          var % dimids(3)       = output_file % dimid_time
     
    23772847          var % task            = "interpolate_2d"
    23782848
    2379        CASE ( 'left scalar', 'right scalar') ! same as right
    2380           var % lod             = 2
     2849       CASE( 'left scalar', 'right scalar') ! same as right
     2850          var % lod             = -1
    23812851          var % ndim            = 3
    23822852          var % dimids(3)       = output_file % dimid_time
     
    23892859          var % task            = "interpolate_3d"
    23902860
    2391        CASE ( 'north scalar', 'south scalar') ! same as south
    2392           var % lod             = 2
     2861       CASE( 'north scalar', 'south scalar') ! same as south
     2862          var % lod             = -1
    23932863          var % ndim            = 3
    23942864          var % dimids(3)       = output_file % dimid_time
     
    24012871          var % task            = "interpolate_3d"
    24022872
    2403        CASE ( 'top scalar', 'top w' )
    2404           var % lod             = 2
     2873       CASE( 'top scalar', 'top w' )
     2874          var % lod             = -1
    24052875          var % ndim            = 3
    24062876          var % dimids(3)       = output_file % dimid_time
     
    24132883          var % task            = "interpolate_3d"
    24142884
    2415        CASE ( 'left u', 'right u' )
    2416           var % lod             = 2
     2885       CASE( 'left u', 'right u' )
     2886          var % lod             = -1
    24172887          var % ndim            = 3
    24182888          var % dimids(3)       = output_file % dimid_time
     
    24252895          var % task            = "interpolate_3d"
    24262896
    2427        CASE ( 'north u', 'south u' )
    2428           var % lod             = 2
     2897       CASE( 'north u', 'south u' )
     2898          var % lod             = -1
    24292899          var % ndim            = 3
    24302900          var % dimids(3)       = output_file % dimid_time    !t
     
    24372907          var % task            = "interpolate_3d"
    24382908
    2439        CASE ( 'top u' )
    2440           var % lod             = 2
     2909       CASE( 'top u' )
     2910          var % lod             = -1
    24412911          var % ndim            = 3
    24422912          var % dimids(3)       = output_file % dimid_time    !t
     
    24492919          var % task            = "interpolate_3d"
    24502920
    2451        CASE ( 'left v', 'right v' )
    2452           var % lod             = 2
     2921       CASE( 'left v', 'right v' )
     2922          var % lod             = -1
    24532923          var % ndim            = 3
    24542924          var % dimids(3)       = output_file % dimid_time
     
    24612931          var % task            = "interpolate_3d"
    24622932
    2463        CASE ( 'north v', 'south v' )
    2464           var % lod             = 2
     2933       CASE( 'north v', 'south v' )
     2934          var % lod             = -1
    24652935          var % ndim            = 3
    24662936          var % dimids(3)       = output_file % dimid_time    !t
     
    24732943          var % task            = "interpolate_3d"
    24742944
    2475        CASE ( 'top v' )
    2476           var % lod             = 2
     2945       CASE( 'top v' )
     2946          var % lod             = -1
    24772947          var % ndim            = 3
    24782948          var % dimids(3)       = output_file % dimid_time    !t
     
    24852955          var % task            = "interpolate_3d"
    24862956
    2487        CASE ( 'left w', 'right w')
    2488           var % lod             = 2
     2957       CASE( 'left w', 'right w')
     2958          var % lod             = -1
    24892959          var % ndim            = 3
    24902960          var % dimids(3)       = output_file % dimid_time
     
    24972967          var % task            = "interpolate_3d"
    24982968
    2499        CASE ( 'north w', 'south w' )
    2500           var % lod             = 2
     2969       CASE( 'north w', 'south w' )
     2970          var % lod             = -1
    25012971          var % ndim            = 3
    25022972          var % dimids(3)       = output_file % dimid_time    !t
     
    25092979          var % task            = "interpolate_3d"
    25102980
    2511        CASE ( 'time series' )
     2981       CASE( 'time series' )
    25122982          var % lod             = 0
    25132983          var % ndim            = 1
     
    25172987          var % task            = "average scalar"
    25182988
    2519        CASE ( 'profile' )
    2520           var % lod             = 2
     2989       CASE( 'constant scalar profile' )
     2990          var % lod             = -1
    25212991          var % ndim            = 2
    25222992          var % dimids(2)       = output_file % dimid_time    !t
     
    25252995          var % dimvarids(1)    = output_file % dimvarids_scl(3)
    25262996          var % to_be_processed = .TRUE.
    2527           var % task            = "profile"
     2997          var % task            = "set profile"
     2998
     2999       CASE( 'large-scale scalar forcing' )
     3000          var % lod             = -1
     3001          var % ndim            = 2
     3002          var % dimids(2)       = output_file % dimid_time    !t
     3003          var % dimids(1)       = output_file % dimids_scl(3) !z
     3004          var % dimvarids(2)    = output_file % dimvarid_time
     3005          var % dimvarids(1)    = output_file % dimvarids_scl(3)
     3006          var % to_be_processed = ls_forcing_variables_required
     3007          var % task            = "average large-scale profile"
     3008
     3009       CASE( 'large-scale w forcing' )
     3010          var % lod             = -1
     3011          var % ndim            = 2
     3012          var % dimids(2)       = output_file % dimid_time    !t
     3013          var % dimids(1)       = output_file % dimids_vel(3) !z
     3014          var % dimvarids(2)    = output_file % dimvarid_time
     3015          var % dimvarids(1)    = output_file % dimvarids_vel(3)
     3016          var % to_be_processed = ls_forcing_variables_required
     3017          var % task            = "average large-scale profile"
    25283018
    25293019       CASE DEFAULT
     
    25603050
    25613051
    2562     SUBROUTINE input_file_list(start_date_string, start_hour, end_hour,        &
    2563                                step_hour, path, prefix, suffix, file_list)
     3052    SUBROUTINE get_input_file_list(start_date_string, start_hour, end_hour,        &
     3053                                   step_hour, path, prefix, suffix, file_list)
    25643054
    25653055       CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
     
    25683058       CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
    25693059
    2570        INTEGER             ::  number_of_files, hour, i
     3060       INTEGER             ::  number_of_intervals, hour, i
    25713061       CHARACTER(LEN=DATE) ::  date_string
    25723062
    2573        number_of_files = end_hour - start_hour + 1
    2574 
    2575        ALLOCATE( file_list(number_of_files) )
    2576 
    2577        DO i = 1, number_of_files
    2578           hour = start_hour + (i-1) * step_hour
     3063       number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
     3064       ALLOCATE( file_list(number_of_intervals + 1) )
     3065
     3066       DO i = 0, number_of_intervals
     3067          hour = start_hour + i * step_hour
    25793068          date_string = add_hours_to(start_date_string, hour)
    25803069
    2581           file_list(i) = TRIM(path) // TRIM(prefix) // TRIM(date_string) //    &
    2582                          TRIM(suffix) // '.nc'
    2583           message = "Set up input file name '" // TRIM(file_list(i)) //"'"
     3070          file_list(i+1) = TRIM(path) // TRIM(prefix) // TRIM(date_string) //    &
     3071                           TRIM(suffix) // '.nc'
     3072          message = "Set up input file name '" // TRIM(file_list(i+1)) //"'"
    25843073          CALL report('input_file_list', message)
    25853074       END DO
    25863075
    2587     END SUBROUTINE input_file_list
     3076    END SUBROUTINE get_input_file_list
    25883077
    25893078
     
    26073096       
    26083097       REAL(dp), ALLOCATABLE                       ::  basic_state_pressure(:)
    2609        TYPE(container), ALLOCATABLE                ::  compute_buffer(:)
     3098       TYPE(container), ALLOCATABLE                ::  preprocess_buffer(:)
    26103099       INTEGER                                     ::  hour, dt
    26113100       INTEGER                                     ::  i, j, k
     
    26173106       
    26183107       CASE( 'velocities' )
    2619           ! Allocate a compute puffer with the same number of arrays as the input
    2620           ALLOCATE( compute_buffer( SIZE(input_buffer) ) )
     3108          ! Allocate a compute buffer with the same number of arrays as the input
     3109          ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) )
    26213110
    26223111          ! Allocate u and v arrays with scalar dimensions
     
    26243113          ny = SIZE(input_buffer(1) % array, 2)
    26253114          nz = SIZE(input_buffer(1) % array, 3)
    2626           ALLOCATE( compute_buffer(1) % array(nx, ny, nz) ) ! u buffer
    2627           ALLOCATE( compute_buffer(2) % array(nx, ny, nz) ) ! v buffer
     3115          ALLOCATE( preprocess_buffer(1) % array(nx, ny, nz) ) ! u buffer
     3116          ALLOCATE( preprocess_buffer(2) % array(nx, ny, nz) ) ! v buffer
    26283117
    26293118 CALL run_control('time', 'alloc')
     
    26323121          CALL centre_velocities( u_face = input_buffer(1) % array,            &
    26333122                                  v_face = input_buffer(2) % array,            &
    2634                                   u_centre = compute_buffer(1) % array,        &
    2635                                   v_centre = compute_buffer(2) % array )
     3123                                  u_centre = preprocess_buffer(1) % array,     &
     3124                                  v_centre = preprocess_buffer(2) % array )
    26363125         
    2637           ! rotate U and V to PALM-4U orientation and overwrite U and V with
    2638           ! rotated velocities
    2639           DO k = 1, nz
    2640           DO j = 2, ny
    2641           DO i = 2, nx
    2642              CALL uv2uvrot( urot = compute_buffer(1) % array(i,j,k),           &
    2643                             vrot = compute_buffer(2) % array(i,j,k),           &
    2644                             rlat = cosmo_grid % lat(j-1),                      &
    2645                             rlon = cosmo_grid % lon(i-1),                      &
    2646                             pollat = phi_cn,                                   &
    2647                             pollon = lambda_cn,                                &
    2648                             u = input_buffer(1) % array(i,j,k),                &
    2649                             v = input_buffer(2) % array(i,j,k) )
    2650           END DO
    2651           END DO
    2652           END DO
     3126          cfg % rotation_method = 'rotated-pole'
     3127          SELECT CASE(cfg % rotation_method)
     3128
     3129          CASE('rotated-pole')
     3130             ! rotate U and V to PALM-4U orientation and overwrite U and V with
     3131             ! rotated velocities
     3132             DO k = 1, nz
     3133             DO j = 2, ny
     3134             DO i = 2, nx
     3135                CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k),     &
     3136                               vrot = preprocess_buffer(2) % array(i,j,k),     &
     3137                               rlat = cosmo_grid % lat(j-1),                   &
     3138                               rlon = cosmo_grid % lon(i-1),                   &
     3139                               pollat = phi_cn,                                &
     3140                               pollon = lambda_cn,                             &
     3141                               u = input_buffer(1) % array(i,j,k),             &
     3142                               v = input_buffer(2) % array(i,j,k) )
     3143             END DO
     3144             END DO
     3145             END DO
     3146
     3147          CASE DEFAULT
     3148             message = "Rotation method '" // TRIM(cfg % rotation_method) //   &
     3149                "' not recognized."
     3150             CALL abort('preprocess', message)
     3151
     3152          END SELECT
     3153
     3154          ! set values
    26533155          input_buffer(1) % array(1,:,:) = 0.0_dp
    26543156          input_buffer(2) % array(1,:,:) = 0.0_dp
     
    26593161 CALL run_control('time', 'comp')
    26603162
    2661           DEALLOCATE( compute_buffer )
     3163          DEALLOCATE( preprocess_buffer )
    26623164 CALL run_control('time', 'alloc')
    26633165
  • palm/trunk/UTIL/inifor/src/inifor.f90

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
     23! Introduced new PALM grid stretching
     24! Renamend initial-condition mode variable 'mode' to 'ic_mode'
     25! Improved log messages
    2326!
    2427!
     
    4851        ONLY:  setup_parameters, setup_grids, setup_variable_tables,           &
    4952               setup_io_groups, fini_grids, fini_variables, fini_io_groups,    &
    50                fini_file_lists, preprocess,                                    &
     53               fini_file_lists, preprocess, origin_lon, origin_lat,            &
    5154               output_file, io_group_list, output_var_table,                   &
    52                cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, mode,            &
    53                imin, imax, jmin, jmax
     55               cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, cfg,             &
     56               average_imin, average_imax, average_jmin, average_jmax
    5457
    5558    USE io
     
    8487
    8588    ! Initialize the netCDF output file and define dimensions
    86     CALL setup_netcdf_dimensions(output_file, palm_grid)
     89    CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date,    &
     90                                 origin_lon, origin_lat)
    8791 CALL run_control('time', 'write')
    8892
    8993    ! Set up the tables containing the input and output variables and set
    9094    ! the corresponding netCDF dimensions for each output variable
    91     CALL setup_variable_tables(mode)
     95    CALL setup_variable_tables(cfg % ic_mode)
    9296 CALL run_control('time', 'write')
    9397
     
    9599    CALL setup_netcdf_variables(output_file % name, output_var_table)
    96100
    97     CALL setup_io_groups() 
     101    CALL setup_io_groups()
    98102 CALL run_control('time', 'init')
    99103
     
    118122 CALL run_control('time', 'comp')
    119123
     124             !TODO: move this assertion into 'preprocess'.
    120125             IF ( .NOT. ALL(input_buffer(:) % is_preprocessed .AND. .TRUE.) )  THEN
    121126                message = "Input buffers for group '" // TRIM(group % kind) // &
     
    159164                      CASE DEFAULT
    160165
    161                           CALL abort("main loop", 'Not a soil variable')
     166                          message = "'" // TRIM(output_var % kind) // "' is not a soil variable"
     167                          CALL abort("main loop", message)
    162168
    163169                      END SELECT
     
    173179                      ALLOCATE( output_arr( 0:output_var % grid % nx,          &
    174180                                            0:output_var % grid % ny,          &
    175                                             0:output_var % grid % nz ) )
     181                                            1:output_var % grid % nz ) )
    176182
    177183 CALL run_control('time', 'alloc')
     
    187193                      ALLOCATE( output_arr( 0:output_var % grid % nx,          &
    188194                                            0:output_var % grid % ny,          &
    189                                             0:output_var % grid % nz ) )
     195                                            1:output_var % grid % nz ) )
    190196 CALL run_control('time', 'alloc')
    191197                     
     
    193199                      CALL average_profile(                                    &
    194200                         input_buffer(output_var % input_id) % array(:,:,:),   &
    195                          output_arr(:,:,:), imin, imax, jmin, jmax,            &
     201                         output_arr(:,:,:), average_imin, average_imax,        &
     202                         average_jmin, average_jmax,                           &
    196203                         output_var % intermediate_grid,                       &
    197204                         output_var % grid)
     
    205212 CALL run_control('time', 'comp')
    206213
    207                    CASE ( 'profile' )
     214                   CASE ( 'set profile' )
    208215                     
    209                       ALLOCATE( output_arr( 1, 1, 0:nz ) )
     216                      ALLOCATE( output_arr( 1, 1, 1:nz ) )
    210217 CALL run_control('time', 'alloc')
    211218
     
    217224                      CASE('ls_forcing_vg')
    218225                          output_arr(1, 1, :) = vg
     226
     227                      CASE('nudging_tau')
     228                          output_arr(1, 1, :) = NUDGING_TAU
    219229
    220230                      CASE DEFAULT
     
    225235                      END SELECT
    226236 CALL run_control('time', 'comp')
     237
     238                   CASE('average large-scale profile')
     239                      message = "Averaging of large-scale forcing profiles " //&
     240                                "has not been implemented, yet."
     241                      CALL abort('main loop', message)
     242                      !ALLOCATE( output_arr( 1, 1, 1:nz ) )
    227243
    228244                   CASE DEFAULT
     
    269285       ELSE
    270286
    271           message = "Skipping IO group '" // TRIM(group % kind) // "'"
     287          message = "Skipping IO group " // TRIM(str(igroup)) // " '" // TRIM(group % kind) // "'"
    272288          IF ( ALLOCATED(group % in_var_list) )  THEN
    273289              message = TRIM(message) // " with input variable '" //           &
     
    291307 CALL run_control('report', 'void')
    292308
    293     message = "Finished writing forcing file '" // TRIM(output_file % name) // &
     309    message = "Finished writing dynamic driver '" // TRIM(output_file % name) // &
    294310              "' successfully."
    295311    CALL report('main loop', message)
  • palm/trunk/UTIL/inifor/src/io.f90

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
     23! Introduced new PALM grid stretching
     24! Updated variable names and metadata for PIDS v1.9 compatibility
     25! Improved handling of the start date string
     26! Better compatibility with older Intel compilers:
     27! - avoiding implicit array allocation with new get_netcdf_variable()
     28!   subroutine instead of function
     29! Improved command line interface:
     30! - Added configuration validation
     31! - New options to configure input file prefixes
     32! - GNU-style short and long option names
     33! - Added version and copyright output
    2334!
    24 ! 
     35!
    2536! Former revisions:
    2637! -----------------
     
    4354    USE control
    4455    USE defs,                                                                  &
    45         ONLY:  DATE, SNAME, PATH, PI, dp, TO_RADIANS, TO_DEGREES, VERSION
     56        ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
    4657    USE netcdf
    4758    USE types
    4859    USE util,                                                                  &
    49         ONLY:  reverse, str
     60        ONLY:  reverse, str, real_to_str
    5061
    5162    IMPLICIT NONE
    5263
     64    INTERFACE get_netcdf_variable
     65        MODULE PROCEDURE get_netcdf_variable_int
     66        MODULE PROCEDURE get_netcdf_variable_real
     67    END INTERFACE get_netcdf_variable
     68
     69    PRIVATE ::  get_netcdf_variable_int, get_netcdf_variable_real
     70
    5371 CONTAINS
     72
     73    SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
     74
     75       CHARACTER(LEN=PATH), INTENT(IN)         ::  in_file
     76       TYPE(nc_var), INTENT(INOUT)             ::  in_var
     77       INTEGER(hp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
     78
     79       INCLUDE 'get_netcdf_variable.inc'
     80
     81    END SUBROUTINE get_netcdf_variable_int
     82
     83
     84    SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
     85
     86       CHARACTER(LEN=PATH), INTENT(IN)      ::  in_file
     87       TYPE(nc_var), INTENT(INOUT)          ::  in_var
     88       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
     89
     90       INCLUDE 'get_netcdf_variable.inc'
     91
     92    END SUBROUTINE get_netcdf_variable_real
     93
    5494
    5595    SUBROUTINE netcdf_define_variable(var, ncid)
     
    5999
    60100        CALL check(nf90_def_var(ncid, var % name, NF90_FLOAT,       var % dimids(1:var % ndim), var % varid))
    61         CALL check(nf90_put_att(ncid, var % varid, "standard_name", var % standard_name))
    62101        CALL check(nf90_put_att(ncid, var % varid, "long_name",     var % long_name))
    63102        CALL check(nf90_put_att(ncid, var % varid, "units",         var % units))
    64         CALL check(nf90_put_att(ncid, var % varid, "lod",           var % lod))
     103        IF ( var % lod .GE. 0 )  THEN
     104           CALL check(nf90_put_att(ncid, var % varid, "lod",           var % lod))
     105        END IF
    65106        CALL check(nf90_put_att(ncid, var % varid, "source",        var % source))
    66107        CALL check(nf90_put_att(ncid, var % varid, "_FillValue",    NF90_FILL_REAL))
     
    94135!> parameters for the PALM-4U computational grid.
    95136!------------------------------------------------------------------------------!
    96     SUBROUTINE parse_command_line_arguments( start_date, hhl_file,             &
    97        soiltyp_file, static_driver_file, input_path, output_file,              &
    98        namelist_file, ug, vg, p0, z0, mode )
    99 
    100        CHARACTER(LEN=PATH), INTENT(INOUT)  ::  hhl_file, soiltyp_file,         &
    101            static_driver_file, input_path, output_file, namelist_file
    102        CHARACTER(LEN=SNAME), INTENT(INOUT) ::  mode
    103        REAL(dp), INTENT(INOUT)             ::  ug, vg, p0, z0
    104        CHARACTER(LEN=DATE), INTENT(INOUT)  ::  start_date
    105 
    106        CHARACTER(LEN=PATH)     ::  option, arg
    107        INTEGER                 ::  arg_count, i
     137    SUBROUTINE parse_command_line_arguments( cfg )
     138
     139       TYPE(inifor_config), INTENT(INOUT) ::  cfg
     140
     141       CHARACTER(LEN=PATH)                ::  option, arg
     142       INTEGER                            ::  arg_count, i
    108143
    109144       arg_count = COMMAND_ARGUMENT_COUNT()
     
    111146
    112147          ! Every option should have an argument.
    113           IF ( MOD(arg_count, 2) .NE. 0 )  THEN
    114              message = "Syntax error in command line."
    115              CALL abort('parse_command_line_arguments', message)
    116           END IF
     148          !IF ( MOD(arg_count, 2) .NE. 0 )  THEN
     149          !   message = "Syntax error in command line."
     150          !   CALL abort('parse_command_line_arguments', message)
     151          !END IF
    117152         
    118153          message = "The -clon and -clat command line options are depricated. " // &
    119154             "Please remove them form your inifor command and specify the " // &
    120155             "location of the PALM-4U origin either" // NEW_LINE(' ') // &
    121              "   - by setting the namelist parameters 'origin_lon' and 'origin_lat, or'" // NEW_LINE(' ') // &
     156             "   - by setting the namelist parameters 'longitude' and 'latitude', or" // NEW_LINE(' ') // &
    122157             "   - by providing a static driver netCDF file via the -static command-line option."
    123158
    124           ! Loop through option/argument pairs.
    125           DO i = 1, arg_count, 2
     159          i = 1
     160          DO WHILE (i .LE. arg_count)
    126161
    127162             CALL GET_COMMAND_ARGUMENT( i, option )
    128              CALL GET_COMMAND_ARGUMENT( i+1, arg )
    129163
    130164             SELECT CASE( TRIM(option) )
    131165
    132              CASE( '-date' )
    133                 start_date = TRIM(arg)
     166             CASE( '-date', '-d', '--date' )
     167                CALL get_option_argument( i, arg )
     168                cfg % start_date = TRIM(arg)
    134169
    135170             ! Elevation of the PALM-4U domain above sea level
    136              CASE( '-z0' )
    137                 READ(arg, *) z0
     171             CASE( '-z0', '-z', '--elevation' )
     172                CALL get_option_argument( i, arg )
     173                READ(arg, *) cfg % z0
    138174
    139175             ! surface pressure, at z0
    140              CASE( '-p0' )
    141                 READ(arg, *) p0
    142 
    143              ! surface pressure, at z0
    144              CASE( '-ug' )
    145                 READ(arg, *) ug
    146 
    147              ! surface pressure, at z0
    148              CASE( '-vg' )
    149                 READ(arg, *) vg
    150 
    151              ! Domain centre geographical longitude
    152              CASE( '-clon' )
     176             CASE( '-p0', '-r', '--surface-pressure' )
     177                CALL get_option_argument( i, arg )
     178                READ(arg, *) cfg % p0
     179
     180             ! geostrophic wind in x direction
     181             CASE( '-ug', '-u', '--geostrophic-u' )
     182                CALL get_option_argument( i, arg )
     183                READ(arg, *) cfg % ug
     184
     185             ! geostrophic wind in y direction
     186             CASE( '-vg', '-v', '--geostrophic-v' )
     187                CALL get_option_argument( i, arg )
     188                READ(arg, *) cfg % vg
     189
     190             ! domain centre geographical longitude and latitude
     191             CASE( '-clon', '-clat' )
    153192                CALL abort('parse_command_line_arguments', message)         
    154193                !READ(arg, *) lambda_cg
    155194                !lambda_cg = lambda_cg * TO_RADIANS
    156 
    157              ! Domain centre geographical latitude
    158              CASE( '-clat' )
    159                 CALL abort('parse_command_line_arguments', message)         
    160195                !READ(arg, *) phi_cg
    161196                !phi_cg = phi_cg * TO_RADIANS
    162197
    163              CASE( '-path' )
    164                  input_path = TRIM(arg)
    165 
    166              CASE( '-hhl' )
    167                  hhl_file = TRIM(arg)
    168 
    169              CASE( '-static' )
    170                  static_driver_file = TRIM(arg)
    171 
    172              CASE( '-soil' )
    173                  soiltyp_file = TRIM(arg)
    174 
    175              CASE( '-o' )
    176                 output_file = TRIM(arg)
    177 
    178              CASE( '-n' )
    179                 namelist_file = TRIM(arg)
    180 
    181              ! Initialization mode: 'profile' / 'volume'
    182              CASE( '-mode' )
    183                 mode = TRIM(arg)
    184 
    185                 SELECT CASE( TRIM(mode) )
    186 
    187                 CASE( 'profile' )
    188 
    189                 CASE DEFAULT
    190                    message = "Mode '" // TRIM(mode) // "' is not supported. " //&
    191                              "Currently, '-mode profile' is the only supported option. " //&
    192                              "Select this one or omit the -mode option entirely."
    193                    CALL abort( 'parse_command_line_arguments', message )
    194                 END SELECT
     198             CASE( '-path', '-p', '--path' )
     199                CALL get_option_argument( i, arg )
     200                 cfg % input_path = TRIM(arg)
     201
     202             CASE( '-hhl', '-l', '--hhl-file' )
     203                CALL get_option_argument( i, arg )
     204                 cfg % hhl_file = TRIM(arg)
     205
     206             CASE( '-static', '-t', '--static-driver' )
     207                CALL get_option_argument( i, arg )
     208                 cfg % static_driver_file = TRIM(arg)
     209
     210             CASE( '-soil', '-s', '--soil-file')
     211                CALL get_option_argument( i, arg )
     212                 cfg % soiltyp_file = TRIM(arg)
     213
     214             CASE( '--flow-prefix')
     215                CALL get_option_argument( i, arg )
     216                 cfg % flow_prefix = TRIM(arg)
     217
     218             CASE( '--radiation-prefix')
     219                CALL get_option_argument( i, arg )
     220                 cfg % radiation_prefix = TRIM(arg)
     221
     222             CASE( '--soil-prefix')
     223                CALL get_option_argument( i, arg )
     224                 cfg % soil_prefix = TRIM(arg)
     225
     226             CASE( '--soilmoisture-prefix')
     227                CALL get_option_argument( i, arg )
     228                 cfg % soilmoisture_prefix = TRIM(arg)
     229
     230             CASE( '-o', '--output' )
     231                CALL get_option_argument( i, arg )
     232                cfg % output_file = TRIM(arg)
     233
     234             CASE( '-n', '--namelist' )
     235                CALL get_option_argument( i, arg )
     236                cfg % namelist_file = TRIM(arg)
     237
     238             ! initial condition mode: 'profile' / 'volume'
     239             CASE( '-mode', '-i', '--init-mode' )
     240                CALL get_option_argument( i, arg )
     241                cfg % ic_mode = TRIM(arg)
     242
     243             ! boundary conditions / forcing mode: 'ideal' / 'real'
     244             CASE( '-f', '--forcing-mode' )
     245                CALL get_option_argument( i, arg )
     246                cfg % bc_mode = TRIM(arg)
     247
     248             CASE( '--version' )
     249                CALL print_version()
     250                STOP
     251
     252             CASE( '--help' )
     253                CALL print_version()
     254                PRINT *, ""
     255                PRINT *, "For a list of command-line options have a look at the README file."
     256                STOP
    195257
    196258             CASE DEFAULT
    197                 message = "unknown option '" // TRIM(option(2:)) // "'."
     259                message = "unknown option '" // TRIM(option) // "'."
    198260                CALL abort('parse_command_line_arguments', message)
    199261
    200262             END SELECT
     263
     264             i = i + 1
    201265
    202266          END DO
     
    210274
    211275   END SUBROUTINE parse_command_line_arguments
     276
     277   
     278   SUBROUTINE get_option_argument(i, arg)
     279      CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
     280      INTEGER, INTENT(INOUT)             ::  i
     281
     282      i = i + 1
     283      CALL GET_COMMAND_ARGUMENT(i, arg)
     284
     285   END SUBROUTINE
     286
     287
     288   SUBROUTINE validate_config(cfg)
     289      TYPE(inifor_config), INTENT(IN) ::  cfg
     290      LOGICAL                         ::  all_files_present
     291
     292      all_files_present = .TRUE.
     293      all_files_present = all_files_present .AND. file_present(cfg % hhl_file)
     294      all_files_present = all_files_present .AND. file_present(cfg % namelist_file)
     295      all_files_present = all_files_present .AND. file_present(cfg % output_file)
     296      all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file)
     297
     298      ! Only check optional static driver file name, if it has been given.
     299      IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
     300         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file)
     301      END IF
     302
     303      IF (.NOT. all_files_present)  THEN
     304         message = "INIFOR configuration invalid; some input files are missing."
     305         CALL abort( 'validate_config', message )
     306      END IF
     307     
     308     
     309      SELECT CASE( TRIM(cfg % ic_mode) )
     310      CASE( 'profile', 'volume')
     311      CASE DEFAULT
     312         message = "Initialization mode '" // TRIM(cfg % ic_mode) //&
     313                   "' is not supported. " //&
     314                   "Please select either 'profile' or 'volume', " //&
     315                   "or omit the -i/--init-mode/-mode option entirely, which corresponds "//&
     316                   "to the latter."
     317         CALL abort( 'validate_config', message )
     318      END SELECT
     319
     320
     321      SELECT CASE( TRIM(cfg % bc_mode) )
     322      CASE( 'real', 'ideal')
     323      CASE DEFAULT
     324         message = "Forcing mode '" // TRIM(cfg % bc_mode) //&
     325                   "' is not supported. " //&
     326                   "Please select either 'real' or 'ideal', " //&
     327                   "or omit the -f/--forcing-mode option entirely, which corresponds "//&
     328                   "to the latter."
     329         CALL abort( 'validate_config', message )
     330      END SELECT
     331
     332
     333   END SUBROUTINE validate_config
     334
     335
     336   LOGICAL FUNCTION file_present(filename)
     337      CHARACTER(LEN=PATH), INTENT(IN) ::  filename
     338
     339      INQUIRE(FILE=filename, EXIST=file_present)
     340
     341      IF (.NOT. file_present)  THEN
     342         message = "The given file '" // "' does not exist."
     343         CALL report('file_present', message)
     344      END IF
     345
     346   END FUNCTION file_present
    212347
    213348
     
    222357!> writes the actual data.
    223358!------------------------------------------------------------------------------!
    224     SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid)
     359   SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid,                  &
     360                                      start_date_string, origin_lon, origin_lat)
    225361
    226362       TYPE(nc_file), INTENT(INOUT)      ::  output_file
    227363       TYPE(grid_definition), INTENT(IN) ::  palm_grid
    228 
    229        CHARACTER (LEN=SNAME) ::  date
     364       CHARACTER (LEN=DATE), INTENT(IN)  ::  start_date_string
     365       REAL(dp), INTENT(IN)              ::  origin_lon, origin_lat
     366
     367       CHARACTER (LEN=8)     ::  date_string
     368       CHARACTER (LEN=10)    ::  time_string
     369       CHARACTER (LEN=5)     ::  zone_string
     370       CHARACTER (LEN=SNAME) ::  history_string
    230371       INTEGER               ::  ncid, nx, ny, nz, nt, dimids(3), dimvarids(3)
    231372       REAL(dp)              ::  z0
    232373
     374       message = "Initializing PALM-4U dynamic driver file '" //               &
     375                 TRIM(output_file % name) // "' and setting up dimensions."
     376       CALL report('setup_netcdf_dimensions', message)
     377
    233378       ! Create the NetCDF file. NF90_CLOBBER selects overwrite mode.
     379#if defined( __netcdf4 )
    234380       CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
     381#else
     382       CALL check(nf90_create(TRIM(output_file % name), NF90_CLOBBER, ncid))
     383#endif
    235384
    236385!
     
    238387!- Section 1: Write global NetCDF attributes
    239388!------------------------------------------------------------------------------
    240        CALL date_and_time(date)
     389       CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
     390       history_string =                                                        &
     391           'Created on '// date_string      //                                 &
     392           ' at '       // time_string(1:2) // ':' // time_string(3:4) //      &
     393           ' (UTC'      // zone_string // ')'
     394
    241395       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title',          'PALM input file for scenario ...'))
    242396       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution',    'Deutscher Wetterdienst, Offenbach'))
    243397       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author',         'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
    244        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        'Created on '//date))
     398       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        TRIM(history_string)))
    245399       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references',     '--'))
    246400       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment',        '--'))
    247        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     '--'))
    248        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     '--'))
     401       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
     402       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
    249403       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
    250404       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
     
    267421          CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
    268422          CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
    269           CALL check( nf90_def_dim(ncid, "z", nz+1, dimids(3)) )
     423          CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
    270424          output_file % dimids_scl = dimids ! save dimids for later
    271425
     
    285439
    286440       ! overwrite third dimid with the one of depth
    287        CALL check(nf90_def_dim(ncid, "depth", SIZE(palm_grid % depths), dimids(3)) )
     441       CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
    288442       output_file % dimids_soil = dimids ! save dimids for later
    289443
    290444       ! overwrite third dimvarid with the one of depth
    291        CALL check(nf90_def_var(ncid, "depth", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
     445       CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
    292446       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
    293447       CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
     
    301455          CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
    302456          CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
    303           CALL check(nf90_def_dim(ncid, "zw", nz, dimids(3)) )
     457          CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
    304458       output_file % dimids_vel = dimids ! save dimids for later
    305459
     
    328482       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "standard_name", "time"))
    329483       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "long_name", "time"))
    330        CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units", "seconds since..."))
     484       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units",     &
     485                               "seconds since " // start_date_string // " UTC"))
    331486
    332487       CALL check(nf90_enddef(ncid))
     
    363518       INTEGER                              ::  i, ncid
    364519
    365        message = "Initializing PALM-4U forcing file '" // TRIM(filename) // "'."
     520       message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
    366521       CALL report('setup_netcdf_variables', message)
    367522
     
    374529
    375530          IF ( var % to_be_processed )  THEN
    376              message = "Defining variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
     531             message = "  variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
    377532             CALL report('setup_netcdf_variables', message)
    378533
     
    386541       CALL check(nf90_close(ncid))
    387542
    388        message = "Forcing file '" // TRIM(filename) // "' initialized successfully."
     543       message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
    389544       CALL report('setup_netcdf_variables', message)
    390545
     
    447602
    448603          input_var => group % in_var_list(1)
    449           buffer(buf_id) % array = get_netcdf_variable( input_file, input_var ) 
     604          CALL get_netcdf_variable(input_file, input_var, buffer(buf_id) % array)
    450605          CALL report('read_input_variables', "Read accumulated " // TRIM(group % in_var_list(1) % name))
    451606
     
    472627             END IF
    473628
    474              buffer(ivar) % array = get_netcdf_variable( input_file, input_var ) 
     629             CALL get_netcdf_variable(input_file, input_var, buffer(ivar) % array)
    475630
    476631             IF ( input_var % is_upside_down )  CALL reverse(buffer(ivar) % array)
     
    545700
    546701          CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
     702          CALL check(nf90_close(ncid))
    547703
    548704       ELSE
     
    555711
    556712    END FUNCTION get_netcdf_attribute
    557 
    558 
    559 
    560     FUNCTION get_netcdf_variable(in_file, in_var) RESULT(buffer)
    561 
    562        CHARACTER(LEN=PATH), INTENT(IN)      ::  in_file
    563        TYPE(nc_var), INTENT(INOUT)          ::  in_var
    564        REAL(dp), ALLOCATABLE                ::  buffer(:,:,:)
    565        INTEGER                              ::  i, ncid, start(3)
    566 
    567 
    568        ! Read in_var NetCDF attributes
    569        IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
    570             nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR )  THEN
    571 
    572           CALL check(nf90_get_att(ncid, in_var % varid, "long_name", in_var % long_name))
    573           CALL check(nf90_get_att(ncid, in_var % varid, "units", in_var % units))
    574 
    575           ! Read in_var NetCDF dimensions
    576           CALL check(nf90_inquire_variable( ncid, in_var % varid,              &
    577                                             ndims  = in_var % ndim,            &
    578                                             dimids = in_var % dimids ))
    579                                        
    580           DO i = 1, in_var % ndim
    581              CALL check(nf90_inquire_dimension( ncid, in_var % dimids(i),      &
    582                                                 name = in_var % dimname(i),    &
    583                                                 len  = in_var % dimlen(i) ))
    584           END DO
    585 
    586           start = (/ 1, 1, 1 /)
    587           IF ( TRIM(in_var % name) .EQ. 'T_SO' )  THEN
    588              ! Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
    589              in_var % dimlen(3) = in_var % dimlen(3) - 1
    590 
    591              ! Start reading from second level, e.g. depth = 0.005 instead of 0.0
    592              start(3) = 2
    593           END IF
    594 
    595           SELECT CASE(in_var % ndim)
    596 
    597           CASE (2)
    598 
    599              ALLOCATE( buffer( in_var % dimlen(1),                             &
    600                                in_var % dimlen(2),                             &
    601                                1 ) )
    602 
    603           CASE (3)
    604 
    605              ALLOCATE( buffer( in_var % dimlen(1),                             &
    606                                in_var % dimlen(2),                             &
    607                                in_var % dimlen(3) ) )
    608           CASE (4)
    609 
    610              ALLOCATE( buffer( in_var % dimlen(1),                             &
    611                                in_var % dimlen(2),                             &
    612                                in_var % dimlen(3) ) )
    613           CASE DEFAULT
    614 
    615              message = "Failed reading NetCDF variable " //                    &
    616                 TRIM(in_var % name) // " with " // TRIM(str(in_var%ndim)) //   &
    617                 " dimensions because only two- and and three-dimensional" //   &
    618                 " variables are supported."
    619              CALL abort('get_netcdf_variable', message)
    620 
    621           END SELECT
    622  CALL run_control('time', 'alloc')
    623          
    624           ! TODO: Check for matching dimensions of buffer and var
    625           CALL check(nf90_get_var( ncid, in_var % varid, buffer,               &
    626                                    start = start,                              &
    627                                    count = in_var % dimlen(1:3) ) )
    628 
    629  CALL run_control('time', 'read')
    630        ELSE
    631 
    632           message = "Failed to read '" // TRIM(in_var % name) // &
    633              "' from file '" // TRIM(in_file) // "'."
    634           CALL report('get_netcdf_variable', message)
    635 
    636        END IF
    637 
    638        CALL check(nf90_close(ncid))
    639 
    640  CALL run_control('time', 'read')
    641 
    642     END FUNCTION get_netcdf_variable
    643713
    644714
     
    657727
    658728       ! Skip time dimension for output
    659        IF ( var_is_time_dependent )  THEN
    660           ndim = var % ndim - 1
    661        ELSE
    662           ndim = var % ndim
    663        END IF
     729       ndim = var % ndim
     730       IF ( var_is_time_dependent )  ndim = var % ndim - 1
    664731
    665732       start(:)      = (/1,1,1,1/)
     
    733800                                   start=start(1:ndim+1) ) )
    734801
    735        CASE ( 'profile' )
     802       CASE ( 'constant scalar profile' )
    736803
    737804          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
    738805                                   start=start(1:ndim+1),                      &
    739806                                   count=count(1:ndim) ) )
     807
     808       CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
     809
     810           message = "Doing nothing in terms of writing large-scale forings."
     811           CALL report('update_output', message)
    740812
    741813       CASE DEFAULT
  • palm/trunk/UTIL/inifor/src/transform.f90

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
     23! Introduced new PALM grid stretching
     24! Removed unnecessary subroutine parameters
     25! Renamed kcur to k_intermediate
    2326!
    2427!
     
    8083       TYPE(grid_definition), INTENT(IN) ::  outgrid
    8184       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
    82        REAL(dp), INTENT(OUT)             ::  out_arr(0:,0:,0:)
    83 
    84        INTEGER :: i, j, k, l, nx, ny, nz
    85 
    86        nx = UBOUND(out_arr, 1)
    87        ny = UBOUND(out_arr, 2)
     85       REAL(dp), INTENT(OUT)             ::  out_arr(0:,0:,:)
     86
     87       INTEGER :: i, j, k, l, nz
     88
    8889       nz = UBOUND(out_arr, 3)
    8990
    90        DO j = 0, ny
    91        DO i = 0, nx
    92        DO k = nz, 0, -1
     91       DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2)
     92       DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1)
     93       DO k = nz, LBOUND(out_arr, 3), -1
    9394
    9495          ! TODO: Remove IF clause and extrapolate based on a critical vertical
     
    101102             out_arr(i,j,k) = 0.0_dp
    102103             DO l = 1, 2
    103                 out_arr(i,j,k) = out_arr(i,j,k) +                                 &
    104                     outgrid % w_verti(i,j,k,l) *                                  &
     104                out_arr(i,j,k) = out_arr(i,j,k) +                              &
     105                    outgrid % w_verti(i,j,k,l) *                               &
    105106                    in_arr(i,j,outgrid % kk(i,j,k, l) )
    106107             END DO
     
    139140    ! I index 0-based for the indices of the outvar to be consistent with the
    140141    ! outgrid indices and interpolation weights.
    141        TYPE(grid_definition), INTENT(IN) ::  outgrid
    142        REAL(dp), INTENT(IN)              ::  invar(0:,0:,0:)
    143        REAL(dp), INTENT(OUT)             ::  outvar(0:,0:,0:)
     142       TYPE(grid_definition), INTENT(IN)  ::  outgrid
     143       REAL(dp), INTENT(IN)               ::  invar(0:,0:,0:)
     144       REAL(dp), INTENT(OUT)              ::  outvar(0:,0:,0:)
    144145       TYPE(nc_var), INTENT(IN), OPTIONAL ::  ncvar
    145146
     
    413414
    414415    END SUBROUTINE rotate_to_cosmo
     416       
    415417
    416418
     
    427429!>                     -------------
    428430!>           jj, lat
    429 !>              ^      j
    430 !>              |       \          i
     431!>              ^        j
     432!>              |         \          i
    431433!>  jj(i,j,2/3) + ... 2 ---\--------/------ 3
    432434!>              |     | ^   \      /        |
     
    459461!>
    460462!------------------------------------------------------------------------------!
    461     SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, cosmo_dxi,     &
    462        cosmo_dyi, palm_clat, palm_clon, palm_ii, palm_jj)
     463    SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon,                &
     464                                          palm_clat, palm_clon,                &
     465                                          palm_ii, palm_jj)
    463466
    464467       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
    465468       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
    466        REAL(dp), INTENT(IN)                       ::  cosmo_dxi, cosmo_dyi
     469       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
    467470       INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT)  ::  palm_ii, palm_jj
    468471
     
    472475       lon0 = cosmo_lon(0)
    473476       lat0 = cosmo_lat(0)
     477       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
     478       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
    474479
    475480       DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
     
    508513       TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
    509514
    510        INTEGER  ::  i, j, k, nx, ny, nz, nlev, kcur
     515       INTEGER  ::  i, j, k, nx, ny, nz, nlev, k_intermediate
    511516       LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
    512517                    point_is_in_current_cell
     
    523528       DO j = 0, ny
    524529
    525           kcur = 0
     530          k_intermediate = 0
    526531
    527532          column_base = palm_intermediate % h(i,j,0)
     
    532537          ! cell, or above column_top. Keep increasing current cell index until
    533538          ! the current cell overlaps with the current_height.
    534           DO k = 0, nz
     539          DO k = 1, nz
    535540
    536541             ! Memorize the top and bottom boundaries of the coarse cell and the
    537542             ! current height within it
    538543             current_height = palm_grid % z(k) + palm_grid % z0
    539              h_top    = palm_intermediate % h(i,j,kcur+1)
    540              h_bottom = palm_intermediate % h(i,j,kcur)
     544             h_top    = palm_intermediate % h(i,j,k_intermediate+1)
     545             h_bottom = palm_intermediate % h(i,j,k_intermediate)
    541546
    542547             point_is_above_grid = (current_height > column_top) !22000m, very unlikely
     
    556561                palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
    557562
     563                message = "PALM-4U grid extends above COSMO-DE model top."
     564                CALL abort('find_vertical_neighbours_and_weights', message)
     565
    558566             ELSE IF (point_is_below_grid)  THEN
    559567
     
    564572                ! cycle through intermediate levels until current
    565573                ! intermediate-grid cell overlaps with current_height
    566                 DO WHILE (.NOT. point_is_in_current_cell .AND. kcur <= nlev-1)
    567                    kcur = kcur + 1
    568 
    569                    h_top    = palm_intermediate % h(i,j,kcur+1)
    570                    h_bottom = palm_intermediate % h(i,j,kcur)
     574                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
     575                   k_intermediate = k_intermediate + 1
     576
     577                   h_top    = palm_intermediate % h(i,j,k_intermediate+1)
     578                   h_bottom = palm_intermediate % h(i,j,k_intermediate)
    571579                   point_is_in_current_cell = (                                &
    572580                      current_height >= h_bottom .AND.                         &
     
    575583                END DO
    576584
    577                 ! kcur = 48 indicates the last section (indices 48 and 49), i.e.
    578                 ! kcur = 49 is not the beginning of a valid cell.
    579                 IF (kcur > nlev-1)  THEN
    580                    message = "Index " // TRIM(str(kcur)) // " is above intermediate grid range."
     585                ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
     586                ! k_intermediate = 49 is not the beginning of a valid cell.
     587                IF (k_intermediate > nlev-1)  THEN
     588                   message = "Index " // TRIM(str(k_intermediate)) //          &
     589                             " is above intermediate grid range."
    581590                   CALL abort('find_vertical_neighbours', message)
    582591                END IF
    583592   
    584                 palm_grid % kk(i,j,k,1) = kcur
    585                 palm_grid % kk(i,j,k,2) = kcur + 1
     593                palm_grid % kk(i,j,k,1) = k_intermediate
     594                palm_grid % kk(i,j,k,2) = k_intermediate + 1
    586595
    587596                ! copmute vertical weights
     
    643652!         
    644653    SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
    645        cosmo_dxi, cosmo_dyi, palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
     654       palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
    646655       
    647656       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
    648        REAL(dp), INTENT(IN)                       ::  cosmo_dxi, cosmo_dyi
     657       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
    649658       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
    650659       INTEGER, DIMENSION(0:,0:,1:), INTENT(IN)   ::  palm_ii, palm_jj
     
    654663       REAL(dp) ::  wl, wp
    655664       INTEGER  ::  i, j
     665
     666       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
     667       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
    656668
    657669       DO j = 0, UBOUND(palm_clon, 2)
  • palm/trunk/UTIL/inifor/src/types.f90

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
     23! Introduced new PALM grid stretching:
     24! - Converted vertical grid_definition coordinte variables to pointers
     25! Improved command line interface:
     26! - Moved INIFOR configuration into a new derived data type
     27! Removed unnecessary variables
    2328!
    2429!
     
    4146 
    4247 USE defs,                                                                     &
    43     ONLY:  dp, PATH, SNAME, LNAME
     48    ONLY:  dp, DATE, PATH, SNAME, LNAME
    4449 USE netcdf,                                                                   &
    4550    ONLY:  NF90_MAX_VAR_DIMS, NF90_MAX_NAME
    4651
    4752 IMPLICIT NONE
     53
     54 TYPE inifor_config
     55    CHARACTER(LEN=DATE)  ::  start_date           !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation
     56
     57    CHARACTER(LEN=PATH)  ::  input_path           !< Path to the input data file directory
     58    CHARACTER(LEN=PATH)  ::  hhl_file             !< Path to the file containing the COSMO-DE HHL variable (height of half layers, i.e. vertical cell faces)
     59    CHARACTER(LEN=PATH)  ::  namelist_file        !< Path to the PALM-4U namelist file
     60    CHARACTER(LEN=PATH)  ::  output_file          !< Path to the INIFOR output file (i.e. PALM-4U dynamic driver')
     61    CHARACTER(LEN=PATH)  ::  soiltyp_file         !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
     62    CHARACTER(LEN=PATH)  ::  static_driver_file   !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
     63
     64    CHARACTER(LEN=SNAME) ::  flow_prefix          !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
     65    CHARACTER(LEN=SNAME) ::  soil_prefix          !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
     66    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses
     67    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses
     68
     69    CHARACTER(LEN=SNAME) ::  bc_mode
     70    CHARACTER(LEN=SNAME) ::  ic_mode
     71    CHARACTER(LEN=SNAME) ::  rotation_method
     72
     73    REAL(dp)             ::  p0
     74    REAL(dp)             ::  ug
     75    REAL(dp)             ::  vg
     76    REAL(dp)             ::  z0 !< Elevation of the PALM-4U domain above sea level [m]
     77 END TYPE inifor_config
    4878
    4979 TYPE grid_definition
     
    5585    INTEGER, ALLOCATABLE  ::  jj(:,:,:)     !< Given a point (i,j,k) in the PALM-4U grid, jj(i,j,l) gives the y index of the l'th horizontl neighbour on the COSMO-DE grid.
    5686    INTEGER, ALLOCATABLE  ::  kk(:,:,:,:)   !< Given a point (i,j,k) in the PALM-4U grid, kk(i,j,k,l) gives the z index of the l'th vertical neighbour in the intermediate grid.
    57     REAL(dp)              ::  dx            !< grid spacing in the first dimension [m]
    58     REAL(dp)              ::  dy            !< grid spacing in the second dimension [m]
    59     REAL(dp)              ::  dz            !< grid spacing in the third dimension [m]
    60     REAL(dp)              ::  dxi           !< inverse grid spacing in the first dimension [m^-1]
    61     REAL(dp)              ::  dyi           !< inverse grid spacing in the second dimension [m^-1]
    62     REAL(dp)              ::  dzi           !< inverse grid spacing in the third dimension [m^-1]
    6387    REAL(dp)              ::  lx            !< domain length in the first dimension [m]
    6488    REAL(dp)              ::  ly            !< domain length in the second dimension [m]
    65     REAL(dp)              ::  lz            !< domain length in the third dimension [m]
    6689    REAL(dp)              ::  x0            !< x coordinate of PALM-4U domain projection centre, i.e. location of zero distortion
    6790    REAL(dp)              ::  y0            !< y coordinate of PALM-4U domain projection centre, i.e. location of zwro distortion
     
    6992    REAL(dp), ALLOCATABLE ::  x(:)          !< coordinates of cell centers in x direction [m]
    7093    REAL(dp), ALLOCATABLE ::  y(:)          !< coordinates of cell centers in y direction [m]
    71     REAL(dp), ALLOCATABLE ::  z(:)          !< coordinates of cell centers in z direction [m]
     94    REAL(dp), POINTER    ::  z(:)          !< coordinates of cell centers in z direction [m]
    7295    REAL(dp), ALLOCATABLE ::  h(:,:,:)      !< heights grid point for intermediate grids [m]
    7396    REAL(dp), POINTER     ::  hhl(:,:,:)    !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from
     
    7699    REAL(dp), ALLOCATABLE ::  xu(:)         !< coordinates of cell faces in x direction [m]
    77100    REAL(dp), ALLOCATABLE ::  yv(:)         !< coordinates of cell faces in y direction [m]
    78     REAL(dp), ALLOCATABLE ::  zw(:)         !< coordinates of cell faces in z direction [m]
     101    REAL(dp), POINTER    ::  zw(:)         !< coordinates of cell faces in z direction [m]
    79102    REAL(dp), ALLOCATABLE ::  lat(:)        !< rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]
    80103    REAL(dp), ALLOCATABLE ::  lon(:)        !< rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]
     
    89112    REAL(dp), ALLOCATABLE ::  w_horiz(:,:,:)   !< weights for bilinear horizontal interpolation
    90113    REAL(dp), ALLOCATABLE ::  w_verti(:,:,:,:) !< weights for linear vertical interpolation
    91  END TYPE
     114 END TYPE grid_definition
    92115
    93116
     
    103126    INTEGER               ::  dimvarids_soil(3)!< NetCDF IDs of the grid coordinates for soil points x, y, depth
    104127    REAL(dp), POINTER     ::  time(:)       ! vector of output time steps
    105  END TYPE
     128 END TYPE nc_file
    106129
    107130
     
    123146    CHARACTER(LEN=SNAME)                  ::  kind                      !< Kind of grid
    124147    CHARACTER(LEN=SNAME)                  ::  task                      !< Processing task that generates this variable, e.g. 'interpolate_2d' or 'average profile'
    125     LOGICAL                               ::  to_be_processed = .FALSE. !< Inifor flag indicating whether variable shall be processed
    126     LOGICAL                               ::  is_read = .FALSE.         !< Inifor flag indicating whether variable has been read
    127     LOGICAL                               ::  is_upside_down  = .FALSE. !< Inifor flag indicating whether variable shall be processed
     148    LOGICAL                               ::  to_be_processed = .FALSE. !< INIFOR flag indicating whether variable shall be processed
     149    LOGICAL                               ::  is_read = .FALSE.         !< INIFOR flag indicating whether variable has been read
     150    LOGICAL                               ::  is_upside_down  = .FALSE. !< INIFOR flag indicating whether variable shall be processed
    128151    TYPE(grid_definition), POINTER        ::  grid                      !< Pointer to the corresponding output grid
    129152    TYPE(grid_definition), POINTER        ::  intermediate_grid         !< Pointer to the corresponding intermediate grid
  • palm/trunk/UTIL/inifor/src/util.f90

    r2718 r3182  
    2121! Current revisions:
    2222! -----------------
     23! Improved real-to-string conversion
    2324!
    2425!
     
    4344        ONLY :  C_CHAR, C_INT, C_PTR, C_SIZE_T
    4445    USE defs,                                                                  &
    45         ONLY :  dp, PI, DATE
     46        ONLY :  dp, PI, DATE, SNAME
    4647
    4748    IMPLICIT NONE
     
    279280    ! Convert a real number to a string in scientific notation
    280281    ! showing four significant digits.
    281     CHARACTER(LEN=11) FUNCTION real_to_str(val, format)
     282    CHARACTER(LEN=SNAME) FUNCTION real_to_str(val, format)
    282283
    283284        REAL(dp), INTENT(IN)                   ::  val
     
    285286
    286287        IF (PRESENT(format))  THEN
    287            WRITE(real_to_str, TRIM(format)) val
     288           WRITE(real_to_str, format) val
    288289        ELSE
    289290           WRITE(real_to_str, '(E11.4)') val
    290            real_to_str = ADJUSTL(real_to_str)
    291291        END IF
     292        real_to_str = ADJUSTL(real_to_str)
    292293
    293294    END FUNCTION real_to_str
Note: See TracChangeset for help on using the changeset viewer.