Index: palm/trunk/UTIL/inifor/src/control.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/control.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/control.f90	(revision 3182)
@@ -21,4 +21,5 @@
 ! Current revisions:
 ! -----------------
+! Added version and copyright output
 ! 
 ! 
@@ -42,5 +43,5 @@
 
     USE defs,                                                                  &
-        ONLY:  LNAME, dp
+        ONLY:  LNAME, dp, VERSION, COPYRIGHT
 
     USE util,                                                                  &
@@ -95,4 +96,10 @@
 
     END SUBROUTINE abort
+
+
+    SUBROUTINE print_version()
+       PRINT *, "INIFOR " // VERSION
+       PRINT *, COPYRIGHT
+    END SUBROUTINE print_version
 
 
Index: palm/trunk/UTIL/inifor/src/defs.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/defs.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/defs.f90	(revision 3182)
@@ -21,4 +21,8 @@
 ! Current revisions:
 ! -----------------
+! Updated defaults for soil extrapolation steps and nudging time-scale
+! Improved handling of the start date string
+! Added gas constant for water vapor
+! Bumped INIFOR version number
 ! 
 ! 
@@ -62,5 +66,6 @@
  REAL(dp), PARAMETER ::  T_SL = 288.15_dp                     !< Reference temperature for computation of COSMO-DE's basic state pressure [K]
  REAL(dp), PARAMETER ::  BETA = 42.0_dp                       !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic state pressure [K]
- 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]
+ 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]
+ REAL(dp), PARAMETER ::  RV   = 461.51_dp                     !< specific gas constant of water vapor [J/kg/K]
  REAL(dp), PARAMETER ::  G    = 9.80665_dp                    !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s]
  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]
@@ -72,7 +77,10 @@
 
  ! INIFOR parameters
- INTEGER, PARAMETER          ::  FILL_ITERATIONS = 10         !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-]
- REAL(dp), PARAMETER         ::  FORCING_FREQ = 3600.0_dp     !< Reference pressure for potential temperature [Pa]
- CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.1.4'            !< path to script for generating input file names
+ INTEGER, PARAMETER          ::  FILL_ITERATIONS = 5          !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-]
+ INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
+ REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
+ CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.3.0'            !< INIFOR version number
+ CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2018 Leibniz Universitaet Hannover' // &
+     NEW_LINE(' ') // ' Copyright 2017-2018 Deutscher Wetterdienst Offenbach' !< Copyright notice
 
  END MODULE defs
Index: palm/trunk/UTIL/inifor/src/get_netcdf_variable.inc
===================================================================
--- palm/trunk/UTIL/inifor/src/get_netcdf_variable.inc	(revision 3182)
+++ palm/trunk/UTIL/inifor/src/get_netcdf_variable.inc	(revision 3182)
@@ -0,0 +1,115 @@
+!> @file src/get_netcdf_variable.inc
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2017-2018 Leibniz Universitaet Hannover
+! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! Initial revision
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+!
+! 
+!
+! Authors:
+! --------
+! @author Eckhard Kadasch
+!
+! Description:
+! ------------
+!> This file provides the implementation of the io.f90 subroutines
+!> get_netcdf_variable_int() and get_netcdf_variable_real(). Both routines
+!> differ only in the type of the netCDF variable to be read and otherwise share
+!> their functionality. In order to avoid duplicating code, both routines source
+!> the implementation below using Fortran's INCLUDE statement.
+!> 
+!------------------------------------------------------------------------------!
+ INTEGER ::  i, ncid, start(3)
+
+ ! Read in_var NetCDF attributes
+ IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
+      nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR )  THEN
+
+    CALL check(nf90_get_att(ncid, in_var % varid, "long_name", in_var % long_name))
+    CALL check(nf90_get_att(ncid, in_var % varid, "units", in_var % units))
+
+    ! Read in_var NetCDF dimensions
+    CALL check(nf90_inquire_variable( ncid, in_var % varid,              &
+                                      ndims  = in_var % ndim,            &
+                                      dimids = in_var % dimids ))
+                                  
+    DO i = 1, in_var % ndim
+       CALL check(nf90_inquire_dimension( ncid, in_var % dimids(i),      &
+                                          name = in_var % dimname(i),    &
+                                          len  = in_var % dimlen(i) ))
+    END DO
+
+    start = (/ 1, 1, 1 /)
+    IF ( TRIM(in_var % name) .EQ. 'T_SO' )  THEN
+       ! Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
+       in_var % dimlen(3) = in_var % dimlen(3) - 1
+
+       ! Start reading from second level, e.g. depth = 0.005 instead of 0.0
+       start(3) = 2
+    END IF
+
+    SELECT CASE(in_var % ndim)
+
+    CASE (2)
+
+       ALLOCATE( buffer( in_var % dimlen(1),                             &
+                         in_var % dimlen(2),                             &
+                         1 ) )
+
+    CASE (3, 4)
+
+       ALLOCATE( buffer( in_var % dimlen(1),                             &
+                         in_var % dimlen(2),                             &
+                         in_var % dimlen(3) ) )
+
+    CASE DEFAULT
+
+       message = "Failed reading NetCDF variable " //                    &
+          TRIM(in_var % name) // " with " // TRIM(str(in_var%ndim)) //   &
+          " dimensions because only two- and and three-dimensional" //   &
+          " variables are supported."
+       CALL abort('get_netcdf_variable', message) 
+
+    END SELECT
+ CALL run_control('time', 'alloc')
+          
+    ! TODO: Check for matching dimensions of buffer and var
+    CALL check(nf90_get_var( ncid, in_var % varid, buffer,               &
+                             start = start,                              &
+                             count = in_var % dimlen(1:3) ) )
+
+ CALL run_control('time', 'read')
+ ELSE
+
+    message = "Failed to read '" // TRIM(in_var % name) // &
+       "' from file '" // TRIM(in_file) // "'."
+    CALL report('get_netcdf_variable', message)
+
+ END IF
+
+ CALL check(nf90_close(ncid))
+
+ CALL run_control('time', 'read')
Index: palm/trunk/UTIL/inifor/src/grid.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/grid.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/grid.f90	(revision 3182)
@@ -21,5 +21,15 @@
 ! Current revisions:
 ! -----------------
-! 
+! Introduced new PALM grid stretching
+! Updated variable names and metadata for PIDS v1.9 compatibility
+! Better compatibility with older Intel compilers:
+! - avoiding implicit array allocation with new get_netcdf_variable()
+!   subroutine instead of function
+! Improved command line interface:
+! - Produce forcing variables when '-mode profile' is being used
+! - Renamend initial-condition mode variable 'mode' to 'ic_mode'
+! - Improved handling of the start date string
+! Removed unnecessary variables and routines
+!
 ! 
 ! Former revisions:
@@ -47,9 +57,9 @@
     USE defs,                                                                  &
         ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp,     &
-               SNAME, LNAME, PATH, FORCING_FREQ, WATER_ID, FILL_ITERATIONS,    &
+               SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS,    &
                BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L
     USE io,                                                                    &
         ONLY:  get_netcdf_variable, get_netcdf_attribute,                      &
-               parse_command_line_arguments
+               parse_command_line_arguments, validate_config
     USE netcdf,                                                                &
         ONLY:  NF90_MAX_NAME, NF90_MAX_VAR_DIMS
@@ -67,36 +77,44 @@
     SAVE
     
-    REAL(dp) ::  phi_equat = 0.0_dp    !< latitude of rotated equator of COSMO-DE grid [rad]
-    REAL(dp) ::  phi_n     = 0.0_dp    !< latitude of rotated pole of COSMO-DE grid [rad]
-    REAL(dp) ::  lambda_n  = 0.0_dp    !< longitude of rotaded pole of COSMO-DE grid [rad]
-    REAL(dp) ::  phi_c     = 0.0_dp    !< rotated-grid latitude of the center of the PALM domain [rad]
-    REAL(dp) ::  lambda_c  = 0.0_dp    !< rotated-grid longitude of the centre of the PALM domain [rad]
-    REAL(dp) ::  phi_cn    = 0.0_dp    !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
-    REAL(dp) ::  lambda_cn = 0.0_dp    !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
-    REAL(dp) ::  gam       = 0.0_dp    !< angle for working around phirot2phi/rlarot2rla bug
-    REAL(dp) ::  dx        = 0.0_dp    !< PALM-4U grid spacing in x direction [m]
-    REAL(dp) ::  dy        = 0.0_dp    !< PALM-4U grid spacing in y direction [m]
-    REAL(dp) ::  dz        = 0.0_dp    !< PALM-4U grid spacing in z direction [m]
-    REAL(dp) ::  dxi       = 0.0_dp    !< inverse PALM-4U grid spacing in x direction [m^-1]
-    REAL(dp) ::  dyi       = 0.0_dp    !< inverse PALM-4U grid spacing in y direction [m^-1]
-    REAL(dp) ::  dzi       = 0.0_dp    !< inverse PALM-4U grid spacing in z direction [m^-1]
-    REAL(dp) ::  lx        = 0.0_dp    !< PALM-4U domain size in x direction [m]
-    REAL(dp) ::  ly        = 0.0_dp    !< PALM-4U domain size in y direction [m]
-    REAL(dp) ::  lz        = 0.0_dp    !< PALM-4U domain size in z direction [m]
-    REAL(dp) ::  ug        = 0.0_dp    !< geostrophic wind in x direction [m/s]
-    REAL(dp) ::  vg        = 0.0_dp    !< geostrophic wind in y direction [m/s]
-    REAL(dp) ::  p0        = 0.0_dp    !< PALM-4U surface pressure, at z0 [Pa]
-    REAL(dp) ::  x0        = 0.0_dp    !< x coordinate of PALM-4U Earth tangent [m] 
-    REAL(dp) ::  y0        = 0.0_dp    !< y coordinate of PALM-4U Earth tangent [m] 
-    REAL(dp) ::  z0        = 0.0_dp    !< Elevation of the PALM-4U domain above sea level [m]
-    REAL(dp) ::  lonmin    = 0.0_dp    !< Minimunm longitude of COSMO-DE's rotated-pole grid
-    REAL(dp) ::  lonmax    = 0.0_dp    !< Maximum longitude of COSMO-DE's rotated-pole grid
-    REAL(dp) ::  latmin    = 0.0_dp    !< Minimunm latitude of COSMO-DE's rotated-pole grid
-    REAL(dp) ::  latmax    = 0.0_dp    !< Maximum latitude of COSMO-DE's rotated-pole grid
-    REAL(dp) ::  latitude  = 0.0_dp    !< geograpohical latitude of the PALM-4U origin, from inipar namelist [deg]
-    REAL(dp) ::  longitude = 0.0_dp    !< geograpohical longitude of the PALM-4U origin, from inipar namelist [deg]
-    REAL(dp) ::  origin_lat= 0.0_dp    !< geograpohical latitude of the PALM-4U origin, from static driver netCDF file [deg]
-    REAL(dp) ::  origin_lon= 0.0_dp    !< geograpohical longitude of the PALM-4U origin, from static driver netCDF file [deg]
-    REAL(dp) ::  end_time  = 0.0_dp    !< PALM-4U simulation time [s]
+    REAL(dp) ::  phi_equat         = 0.0_dp       !< latitude of rotated equator of COSMO-DE grid [rad]
+    REAL(dp) ::  phi_n             = 0.0_dp       !< latitude of rotated pole of COSMO-DE grid [rad]
+    REAL(dp) ::  lambda_n          = 0.0_dp       !< longitude of rotaded pole of COSMO-DE grid [rad]
+    REAL(dp) ::  phi_c             = 0.0_dp       !< rotated-grid latitude of the center of the PALM domain [rad]
+    REAL(dp) ::  lambda_c          = 0.0_dp       !< rotated-grid longitude of the centre of the PALM domain [rad]
+    REAL(dp) ::  phi_cn            = 0.0_dp       !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
+    REAL(dp) ::  lambda_cn         = 0.0_dp       !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
+    REAL(dp) ::  gam               = 0.0_dp       !< angle for working around phirot2phi/rlarot2rla bug
+    REAL(dp) ::  dx                = 0.0_dp       !< PALM-4U grid spacing in x direction [m]
+    REAL(dp) ::  dy                = 0.0_dp       !< PALM-4U grid spacing in y direction [m]
+    REAL(dp) ::  dz(10)            = -1.0_dp      !< PALM-4U grid spacing in z direction [m]
+    REAL(dp) ::  dz_max            = 1000.0_dp    !< maximum vertical grid spacing [m]
+    REAL(dp) ::  dz_stretch_factor = 1.08_dp      !< factor for vertical grid stretching [m]
+    REAL(dp) ::  dz_stretch_level  = -9999999.9_dp!< height above which the vertical grid will be stretched [m]
+    REAL(dp) ::  dz_stretch_level_start(9) = -9999999.9_dp !< namelist parameter
+    REAL(dp) ::  dz_stretch_level_end(9) = 9999999.9_dp !< namelist parameter
+    REAL(dp) ::  dz_stretch_factor_array(9) = 1.08_dp !< namelist parameter
+    REAL(dp) ::  dxi               = 0.0_dp       !< inverse PALM-4U grid spacing in x direction [m^-1]
+    REAL(dp) ::  dyi               = 0.0_dp       !< inverse PALM-4U grid spacing in y direction [m^-1]
+    REAL(dp) ::  dzi               = 0.0_dp       !< inverse PALM-4U grid spacing in z direction [m^-1]
+    REAL(dp) ::  lx                = 0.0_dp       !< PALM-4U domain size in x direction [m]
+    REAL(dp) ::  ly                = 0.0_dp       !< PALM-4U domain size in y direction [m]
+    REAL(dp) ::  ug                = 0.0_dp       !< geostrophic wind in x direction [m/s]
+    REAL(dp) ::  vg                = 0.0_dp       !< geostrophic wind in y direction [m/s]
+    REAL(dp) ::  p0                = 0.0_dp       !< PALM-4U surface pressure, at z0 [Pa]
+    REAL(dp) ::  x0                = 0.0_dp       !< x coordinate of PALM-4U Earth tangent [m] 
+    REAL(dp) ::  y0                = 0.0_dp       !< y coordinate of PALM-4U Earth tangent [m] 
+    REAL(dp) ::  z0                = 0.0_dp       !< Elevation of the PALM-4U domain above sea level [m]
+    REAL(dp) ::  z_top             = 0.0_dp       !< height of the scalar top boundary [m]
+    REAL(dp) ::  zw_top            = 0.0_dp       !< height of the vertical velocity top boundary [m]
+    REAL(dp) ::  lonmin            = 0.0_dp       !< Minimunm longitude of COSMO-DE's rotated-pole grid
+    REAL(dp) ::  lonmax            = 0.0_dp       !< Maximum longitude of COSMO-DE's rotated-pole grid
+    REAL(dp) ::  latmin            = 0.0_dp       !< Minimunm latitude of COSMO-DE's rotated-pole grid
+    REAL(dp) ::  latmax            = 0.0_dp       !< Maximum latitude of COSMO-DE's rotated-pole grid
+    REAL(dp) ::  latitude          = 0.0_dp       !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
+    REAL(dp) ::  longitude         = 0.0_dp       !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
+    REAL(dp) ::  origin_lat        = 0.0_dp       !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg]
+    REAL(dp) ::  origin_lon        = 0.0_dp       !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg]
+    REAL(dp) ::  rotation_angle    = 0.0_dp       !< clockwise angle the PALM-4U north is rotated away from geographical north [deg]
+    REAL(dp) ::  end_time          = 0.0_dp       !< PALM-4U simulation time [s]
 
     REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hhl             !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
@@ -106,9 +124,19 @@
     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlon            !< longitudes of COSMO-DE's rotated-pole grid
     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlat            !< latitudes of COSMO-DE's rotated-pole grid
-    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  time
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  time            !< output times 
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  x               !< base palm grid x coordinate vector pointed to by grid_definitions
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  xu              !< base palm grid xu coordinate vector pointed to by grid_definitions
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  y               !< base palm grid y coordinate vector pointed to by grid_definitions
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  yv              !< base palm grid yv coordinate vector pointed to by grid_definitions
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z_column        !< base palm grid z coordinate vector including the top boundary coordinate (entire column)
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw_column       !< base palm grid zw coordinate vector including the top boundary coordinate (entire column)
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z               !< base palm grid z coordinate vector pointed to by grid_definitions
+    REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw              !< base palm grid zw coordinate vector pointed to by grid_definitions
 
     INTEGER(hp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  soiltyp      !< COSMO-DE soil type map
+    INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
+    INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
     INTEGER ::  i     !< indexing variable
-    INTEGER ::  imin, imax, jmin, jmax !< index ranges for profile averaging
+    INTEGER ::  average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging
     INTEGER ::  k     !< indexing variable
     INTEGER ::  nt    !< number of output time steps
@@ -130,4 +158,6 @@
     LOGICAL ::  init_variables_required
     LOGICAL ::  boundary_variables_required
+    LOGICAL ::  ls_forcing_variables_required
+    LOGICAL ::  profile_grids_required
 
     INTEGER ::  n_invar = 0 !< number of variables in the input variable table
@@ -193,14 +223,10 @@
     TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
  
-    NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude
+    NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude,             &
+                      dz_max, dz_stretch_factor, dz_stretch_level,             & !< old grid stretching parameters
+                      dz_stretch_level_start, dz_stretch_level_end               !< new grid stretching parameters
     NAMELIST /d3par/  end_time
     
     CHARACTER(LEN=LNAME) ::  nc_source_text     = ''  !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...'
-    CHARACTER(LEN=DATE)  ::  start_date         = ''  !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation
-    CHARACTER(LEN=PATH)  ::  hhl_file           = ''  !< Path to the file containing the COSMO-DE HHL variable (height of half layers, i.e. vertical cell faces)
-    CHARACTER(LEN=PATH)  ::  namelist_file      = ''  !< Path to the PALM-4U namelist file
-    CHARACTER(LEN=PATH)  ::  static_driver_file = ''  !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
-    CHARACTER(LEN=PATH)  ::  soiltyp_file       = ''  !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
-    CHARACTER(LEN=PATH)  ::  input_path         = ''  !< Path to the input data file directory
     CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  flow_files
     CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_moisture_files
@@ -212,7 +238,8 @@
     CHARACTER(LEN=SNAME) ::  radiation_suffix     !< Suffix of radiation input files, e.g. 'radiation'
     CHARACTER(LEN=SNAME) ::  soilmoisture_suffix  !< Suffix of input files for soil moisture spin-up, e.g. 'soilmoisture'
-    CHARACTER(LEN=SNAME) ::  mode                 !< INIFOR's initialization mode, 'profile' or 'volume'
                           
     TYPE(nc_file) ::  output_file
+
+    TYPE(inifor_config) ::  cfg !< Container of INIFOR configuration
 
  CONTAINS
@@ -224,5 +251,5 @@
 ! Section 1: Define default parameters
 !------------------------------------------------------------------------------
-       start_date = '2013072100'
+       cfg % start_date = '2013072100'
        end_hour = 2
        start_hour_soil = -2
@@ -249,10 +276,10 @@
        origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
        origin_lon = 13.082744_dp * TO_RADIANS
-       z0 = 35.0_dp
+       cfg % z0 = 35.0_dp
 
        ! Default atmospheric parameters
        ug = 0.0_dp
        vg = 0.0_dp
-       p0 = P_SL
+       cfg % p0 = P_SL
 
        ! Parameters for file names
@@ -262,6 +289,12 @@
        start_hour_soilmoisture = start_hour_flow - 2
        end_hour_soilmoisture = start_hour_flow
-       step_hour = 1
+       step_hour = FORCING_STEP
+
        input_prefix = 'laf'  ! 'laf' for COSMO-DE analyses
+       cfg % flow_prefix = input_prefix
+       cfg % soil_prefix = input_prefix
+       cfg % radiation_prefix = input_prefix
+       cfg % soilmoisture_prefix  = input_prefix
+
        flow_suffix = '-flow'
        soil_suffix = '-soil'
@@ -273,32 +306,48 @@
 !------------------------------------------------------------------------------
 
-       ! Set default paths
-       input_path         = './'
-       hhl_file           = ''
-       soiltyp_file       = ''
-       namelist_file      = './namelist'
-       output_file % name = './palm-4u-input.nc'
+       ! Set default paths and modes
+       cfg % input_path         = './'
+       cfg % hhl_file           = ''
+       cfg % soiltyp_file       = ''
+       cfg % namelist_file      = './namelist'
+       cfg % static_driver_file = ''
+       cfg % output_file = './palm-4u-input.nc'
+       cfg % ic_mode = 'volume'
+       cfg % bc_mode = 'real'
 
        ! Update default file names and domain centre
-       CALL parse_command_line_arguments( start_date, hhl_file, soiltyp_file,  &
-               static_driver_file, input_path, output_file % name,             &
-               namelist_file, ug, vg, p0, z0, mode )
-
-       init_variables_required     = .TRUE.
-       boundary_variables_required = (TRIM(mode) .NE. 'profile')
-
-       CALL normalize_path(input_path)
-       IF (TRIM(hhl_file) == '')  hhl_file = TRIM(input_path) // 'hhl.nc'
-       IF (TRIM(soiltyp_file) == '')  soiltyp_file = TRIM(input_path) // 'soil.nc'
-
-       CALL report('setup_parameters', "       data path: " // TRIM(input_path))
-       CALL report('setup_parameters', "        hhl file: " // TRIM(hhl_file))
-       CALL report('setup_parameters', "    soiltyp file: " // TRIM(soiltyp_file))
-       CALL report('setup_parameters', "   namelist file: " // TRIM(namelist_file))
-       CALL report('setup_parameters', "output data file: " // TRIM(output_file % name))
+       CALL parse_command_line_arguments( cfg )
+
+       output_file % name = cfg % output_file
+       z0 = cfg % z0
+       p0 = cfg % p0
+
+       init_variables_required = .TRUE.
+       boundary_variables_required = TRIM( cfg % bc_mode ) == 'real'
+       ls_forcing_variables_required = TRIM( cfg % bc_mode ) == 'ideal'
+
+       IF ( ls_forcing_variables_required )  THEN
+          message = "Averaging of large-scale forcing profiles " //            &
+                    "has not been implemented, yet."
+          CALL abort('setup_parameters', message)
+       END IF
+
+       CALL normalize_path(cfg % input_path)
+       IF (TRIM(cfg % hhl_file) == '')  cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
+       IF (TRIM(cfg % soiltyp_file) == '')  cfg % soiltyp_file = TRIM(cfg % input_path) // 'soil.nc'
+
+       CALL validate_config( cfg ) 
+
+       CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode))
+       CALL report('setup_parameters', "       forcing mode: " // TRIM(cfg % bc_mode))
+       CALL report('setup_parameters', "          data path: " // TRIM(cfg % input_path))
+       CALL report('setup_parameters', "           hhl file: " // TRIM(cfg % hhl_file))
+       CALL report('setup_parameters', "       soiltyp file: " // TRIM(cfg % soiltyp_file))
+       CALL report('setup_parameters', "      namelist file: " // TRIM(cfg % namelist_file))
+       CALL report('setup_parameters', "   output data file: " // TRIM(output_file % name))
 
  CALL run_control('time', 'init')
        ! Read in namelist parameters
-       OPEN(10, FILE=namelist_file)
+       OPEN(10, FILE=cfg % namelist_file)
        READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz
        READ(10, NML=d3par)  ! end_time
@@ -306,15 +355,15 @@
  CALL run_control('time', 'read')
 
-       end_hour = CEILING(end_time / FORCING_FREQ)
+       end_hour = CEILING( end_time / 3600.0 * step_hour )
 
        ! Generate input file lists
-       CALL input_file_list(start_date, start_hour_flow, end_hour, step_hour,  &
-          input_path, input_prefix, flow_suffix, flow_files)
-       CALL input_file_list(start_date, start_hour_soil, end_hour, step_hour,  &
-          input_path, input_prefix, soil_suffix, soil_files)
-       CALL input_file_list(start_date, start_hour_radiation, end_hour, step_hour, &
-          input_path, input_prefix, radiation_suffix, radiation_files)
-       CALL input_file_list(start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
-          input_path, input_prefix, soilmoisture_suffix, soil_moisture_files)
+       CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour,  &
+          cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files)
+       CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour,  &
+          cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files)
+       CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, &
+          cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files)
+       CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
+          cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
 
 !
@@ -322,9 +371,5 @@
 ! Section 3: Check for consistency
 !------------------------------------------------------------------------------
-       IF (dx*dy*dz .EQ. 0.0_dp)  THEN
-          message = "Grid cells have zero volume. Grid spacings are probably"//&
-             " set incorrectly in namelist file '" // TRIM(namelist_file) // "'."
-          CALL abort('setup_parameters', message) 
-       END IF
+
 !
 !------------------------------------------------------------------------------
@@ -339,14 +384,14 @@
        ! Read COSMO-DE soil type map
        cosmo_var % name = 'SOILTYP'
-       soiltyp = NINT(get_netcdf_variable(soiltyp_file, cosmo_var), hp)
+       CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp)
 
        message = 'Reading PALM-4U origin from'
-       IF (TRIM(static_driver_file) .NE. '')  THEN
-
-          origin_lon = get_netcdf_attribute(static_driver_file, 'origin_lon')
-          origin_lat = get_netcdf_attribute(static_driver_file, 'origin_lat')
+       IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
+
+          origin_lon = get_netcdf_attribute(cfg % static_driver_file, 'origin_lon')
+          origin_lat = get_netcdf_attribute(cfg % static_driver_file, 'origin_lat')
 
           message = TRIM(message) // " static driver file '"                   &
-                                  // TRIM(static_driver_file) // "'"
+                                  // TRIM(cfg % static_driver_file) // "'"
 
 
@@ -357,5 +402,5 @@
 
           message = TRIM(message) // " namlist file '"                         &
-                                  // TRIM(namelist_file) // "'"
+                                  // TRIM(cfg % namelist_file) // "'"
 
        END IF
@@ -368,5 +413,5 @@
        ! Read in COSMO-DE heights of half layers (vertical cell faces)
        cosmo_var % name = 'HHL'
-       hhl = get_netcdf_variable(hhl_file, cosmo_var)
+       CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl)
  CALL run_control('time', 'read')
 
@@ -392,5 +437,4 @@
        lx = (nx+1) * dx
        ly = (ny+1) * dy
-       lz = (nz+1) * dz
        
        ! PALM-4U point of Earth tangency
@@ -399,5 +443,5 @@
 
        ! time vector
-       nt = CEILING(end_time / FORCING_FREQ) + 1
+       nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1
        ALLOCATE( time(nt) )
        CALL linspace(0.0_dp, 3600.0_dp * (nt-1), time)
@@ -463,13 +507,42 @@
     SUBROUTINE setup_grids() ! setup_grids(inifor_settings(with nx, ny, nz,...))
        CHARACTER ::  interp_mode
-        
+
 !------------------------------------------------------------------------------
-! Section 1: Define model and initialization grids
+! Section 0: Define base PALM-4U grid coordinate vectors
+!------------------------------------------------------------------------------
+       ! palm x y z, we allocate the column to nz+1 in order to include the top
+       ! scalar boundary. The interpolation grids will be associated with
+       ! a shorter column that omits the top element.
+       
+       ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
+       CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x)
+       CALL linspace(0.5_dp * dy, ly - 0.5_dp * dy, y)
+       CALL stretched_z(z_column, dz, dz_max=dz_max,                           &
+                        dz_stretch_factor=dz_stretch_factor,                   &
+                        dz_stretch_level=dz_stretch_level,                     &
+                        dz_stretch_level_start=dz_stretch_level_start,         &
+                        dz_stretch_level_end=dz_stretch_level_end,             &
+                        dz_stretch_factor_array=dz_stretch_factor_array)
+       z(1:nz) = z_column(1:nz)
+       z_top = z_column(nz+1)
+
+       ! palm xu yv zw, compared to the scalar grid, velocity coordinates
+       ! contain one element less.
+       ALLOCATE( xu(1:nx),  yv(1:ny), zw(1:nz-1), zw_column(1:nz))
+       CALL linspace(dx, lx - dx, xu)
+       CALL linspace(dy, ly - dy, yv)
+       CALL midpoints(z_column, zw_column)
+       zw(1:nz-1) = zw_column(1:nz-1)
+       zw_top     = zw_column(nz)
+
+
+!------------------------------------------------------------------------------
+! Section 1: Define initialization and boundary grids
 !------------------------------------------------------------------------------
        CALL init_grid_definition('palm', grid=palm_grid,                       &
                xmin=0.0_dp, xmax=lx,                                           &
                ymin=0.0_dp, ymax=ly,                                           &
-               zmin=0.0_dp, zmax=lz, x0=x0, y0=y0, z0=z0,                      &
-               nx=nx, ny=ny, nz=nz, mode=mode)
+               x0=x0, y0=y0, z0=z0,                      &
+               nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode)
 
        ! Subtracting 1 because arrays will be allocated with nlon + 1 elements.
@@ -477,5 +550,5 @@
                xmin=lonmin, xmax=lonmax,                                       &
                ymin=latmin, ymax=latmax,                                       &
-               zmin=0.0_dp, zmax=51.0_dp, x0=x0, y0=y0, z0=0.0_dp,             &
+               x0=x0, y0=y0, z0=0.0_dp,             &
                nx=nlon-1, ny=nlat-1, nz=nlev-1)
 
@@ -487,5 +560,5 @@
                xmin=0.0_dp, xmax=lx,                                           &
                ymin=0.0_dp, ymax=ly,                                           &
-               zmin=0.0_dp, zmax=lz, x0=x0, y0=y0, z0=z0,                      &
+               x0=x0, y0=y0, z0=z0,                      &
                nx=nx, ny=ny, nz=nlev-2)
 
@@ -493,48 +566,39 @@
                xmin = dx, xmax = lx - dx,                                      &
                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-               zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
                x0=x0, y0=y0, z0 = z0,                                          &
                nx = nx-1, ny = ny, nz = nz,                                    &
-               dx = dx, dy = dy, dz = dz, mode=mode)
+               z=z, ic_mode=cfg % ic_mode)
 
        CALL init_grid_definition('boundary', grid=v_initial_grid,              &
                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                ymin = dy, ymax = ly - dy,                                      &
-               zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
                x0=x0, y0=y0, z0 = z0,                                          &
                nx = nx, ny = ny-1, nz = nz,                                    &
-               dx = dx, dy = dy, dz = dz, mode=mode)
+               z=z, ic_mode=cfg % ic_mode)
 
        CALL init_grid_definition('boundary', grid=w_initial_grid,              &
                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-               zmin = dz, zmax = lz - dz,                                      &
                x0=x0, y0=y0, z0 = z0,                                          &
                nx = nx, ny = ny, nz = nz-1,                                    &
-               dx = dx, dy = dy, dz = dz, mode=mode)
+               z=zw, ic_mode=cfg % ic_mode)
 
        CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate,      &
                xmin = dx, xmax = lx - dx,                                      &
                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-               zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
                x0=x0, y0=y0, z0 = z0,                                          &
-               nx = nx-1, ny = ny, nz = nlev - 2,                              &
-               dx = dx, dy = dy, dz = dz)
+               nx = nx-1, ny = ny, nz = nlev - 2)
 
        CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate,      &
                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                ymin = dy, ymax = ly - dy,                                      &
-               zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
                x0=x0, y0=y0, z0 = z0,                                          &
-               nx = nx, ny = ny-1, nz = nlev - 2,                              &
-               dx = dx, dy = dy, dz = dz)
+               nx = nx, ny = ny-1, nz = nlev - 2)
 
        CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate,      &
                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-               zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
                x0=x0, y0=y0, z0 = z0,                                          &
-               nx = nx, ny = ny, nz = nlev - 2,                                &
-               dx = dx, dy = dy, dz = dz)
+               nx = nx, ny = ny, nz = nlev - 2)
 
       IF (boundary_variables_required)  THEN
@@ -546,320 +610,240 @@
                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.5_dp * dz, zmax = lz - 0.5_dp * dz,                   &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=scalars_west_grid,           &
                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.5_dp * dz, zmax = lz - 0.5_dp * dz,                   &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=scalars_north_grid,          &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=scalars_south_grid,          &
                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
-                  zmin =  0.5_dp * dz, zmax = lz - 0.5_dp * dz,                   &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=scalars_top_grid,            &
                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz,              &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny, nz = 0,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = ny, nz = 1, z=(/z_top/))
 
           CALL init_grid_definition('boundary', grid=u_east_grid,                 &
                   xmin = lx, xmax = lx,                                           &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=u_west_grid,                 &
                   xmin = 0.0_dp, xmax = 0.0_dp,                                   &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=u_north_grid,                &
                   xmin = dx, xmax = lx - dx,                                      &
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx-1, ny = 0, nz = nz,                                     &
-                  dx = dx, dy = dy, dz = dz)
-
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx-1, ny = 0, nz = nz, z=z)
+    
           CALL init_grid_definition('boundary', grid=u_south_grid,                &
                   xmin = dx, xmax = lx - dx,                                      &
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx-1, ny = 0, nz = nz,                                     &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx-1, ny = 0, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=u_top_grid,                  &
                   xmin = dx, xmax = lx - dx,                                      &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin = lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz,               &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx-1, ny = ny, nz = 0,                                     &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx-1, ny = ny, nz = 1, z=(/z_top/))
 
           CALL init_grid_definition('boundary', grid=v_east_grid,                 &
                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
                   ymin = dy, ymax = ly - dy,                                      &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny-1, nz = nz,                                     &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny-1, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=v_west_grid,                 &
                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
                   ymin = dy, ymax = ly - dy,                                      &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny-1, nz = nz,                                     &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny-1, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=v_north_grid,                &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = ly, ymax = ly,                                           &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=v_south_grid,                &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = 0.0_dp, ymax = 0.0_dp,                                   &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                    &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nz,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=v_top_grid,                  &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = dy, ymax = ly - dy,                                      &
-                  zmin = lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz,               &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny-1, nz = 0,                                     &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = ny-1, nz = 1, z=(/z_top/))
 
           CALL init_grid_definition('boundary', grid=w_east_grid,                 &
                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  dz, zmax = lz - dz,                                     &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nz - 1,                                   &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nz - 1, z=zw)
 
           CALL init_grid_definition('boundary', grid=w_west_grid,                 &
                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin = dz, zmax = lz - dz,                                      &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nz - 1,                                   &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nz - 1, z=zw)
 
           CALL init_grid_definition('boundary', grid=w_north_grid,                &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
-                  zmin = dz, zmax = lz - dz,                                      &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nz - 1,                                   &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nz - 1, z=zw)
 
           CALL init_grid_definition('boundary', grid=w_south_grid,                &
                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
-                  zmin = dz, zmax = lz - dz,                                      &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nz - 1,                                   &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nz - 1, z=zw)
 
           CALL init_grid_definition('boundary', grid=w_top_grid,                  &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin = lz, zmax = lz,                                           &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny, nz = 0,                                       &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = ny, nz = 1, z=(/zw_top/))
 
           CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate,   &
                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate,   &
                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate,  &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate,  &
                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate,    &
                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny, nz = nlev - 2,                                &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate,         &
                   xmin = lx, xmax = lx,                                           &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate,         &
                   xmin = 0.0_dp, xmax = 0.0_dp,                                   &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate,        &
                   xmin = dx, xmax = lx - dx,                                      &
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx-1, ny = 0, nz = nlev - 2,                               &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx-1, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate,        &
                   xmin = dx, xmax = lx - dx,                                      &
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx-1, ny = 0, nz = nlev - 2,                               &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx-1, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate,          &
                   xmin = dx, xmax = lx - dx,                                      &
                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx-1, ny = ny, nz = nlev - 2,                              &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx-1, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate,         &
                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
                   ymin = dy, ymax = ly - dy,                                      &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny-1, nz = nlev - 2,                               &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny-1, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate,         &
                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
                   ymin = dy, ymax = ly - dy,                                      &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny-1, nz = nlev - 2,                               &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny-1, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate,        &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = ly, ymax = ly,                                           &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate,        &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = 0.0_dp, ymax = 0.0_dp,                                   &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate,          &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = dy, ymax = ly - dy,                                      &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny-1, nz = nlev - 2,                              &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = ny-1, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate,         &
                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = 0, ny = ny, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = 0, nz = nlev - 2)
 
           CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
-                  zmin =  0.0_dp, zmax = 0.0_dp,                                  &
-                  x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny, nz = nlev - 2,                                &
-                  dx = dx, dy = dy, dz = dz)
+                  x0=x0, y0=y0, z0 = z0,                                          &
+                  nx = nx, ny = ny, nz = nlev - 2)
        END IF
 
@@ -869,36 +853,31 @@
 !------------------------------------------------------------------------------
 
-       IF (TRIM(mode) == 'profile')  THEN
+       profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR.              &
+                                  TRIM(cfg % bc_mode) == 'ideal' )
+
+       IF (profile_grids_required)  THEN
           CALL init_grid_definition('boundary', grid=scalar_profile_grid,      &
                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz,                 &
                   x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nz,                                     &
-                  dx = dx, dy = dy, dz = dz)
+                  nx = 0, ny = 0, nz = nz, z=z)
 
           CALL init_grid_definition('boundary', grid=w_profile_grid,           &
                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  zmin = dz, zmax = lz - dz,                                   &
                   x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nz - 1,                                 &
-                  dx = dx, dy = dy, dz = dz)
+                  nx = 0, ny = 0, nz = nz - 1, z=zw)
 
           CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,&
                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  zmin = 0.0_dp, zmax = 0.0_dp,                                &
                   x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nlev - 2,                               &
-                  dx = dx, dy = dy, dz = dz)
+                  nx = 0, ny = 0, nz = nlev - 2, z=z)
 
           CALL init_grid_definition('boundary', grid=w_profile_intermediate,   &
                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  zmin = 0.0_dp, zmax = 0.0_dp,                                &
                   x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nlev - 2,                               &
-                  dx = dx, dy = dy, dz = dz)
+                  nx = 0, ny = 0, nz = nlev - 2, z=zw)
        END IF
 
@@ -908,5 +887,5 @@
 !------------------------------------------------------------------------------
        interp_mode = 's'
-       CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, mode=mode)
+       CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg % ic_mode)
        IF (boundary_variables_required)  THEN
           CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode)
@@ -918,5 +897,5 @@
 
        interp_mode = 'u'
-       CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, mode=mode)
+       CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
        IF (boundary_variables_required)  THEN
           CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode)
@@ -928,5 +907,5 @@
 
        interp_mode = 'v'
-       CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, mode=mode)
+       CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
        IF (boundary_variables_required)  THEN
           CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode)
@@ -938,5 +917,5 @@
 
        interp_mode = 'w'
-       CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, mode=mode)
+       CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
        IF (boundary_variables_required)  THEN
           CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode)
@@ -947,6 +926,7 @@
        END IF
 
-       IF (TRIM(mode) == 'profile')  THEN
-           CALL setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax) 
+       IF (TRIM(cfg % ic_mode) == 'profile')  THEN
+           CALL setup_averaging(cosmo_grid, palm_intermediate,                 &
+                average_imin, average_imax, average_jmin, average_jmax) 
        END IF
         
@@ -955,10 +935,10 @@
 
 
-    SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, mode)
+    SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode)
 
        TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
        TYPE(grid_definition), INTENT(INOUT), TARGET ::  palm_grid, palm_intermediate
        CHARACTER, INTENT(IN)                        ::  kind
-       CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  mode
+       CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
 
        TYPE(grid_definition), POINTER      ::  grid
@@ -966,10 +946,9 @@
        REAL(dp), DIMENSION(:,:,:), POINTER ::  h
 
-       LOGICAL :: setup_vertical = .TRUE.
-
-       IF (PRESENT(mode))  THEN
-          IF (TRIM(mode) == 'profile')  setup_vertical = .FALSE.
-       ELSE
-          setup_vertical = .TRUE.
+       LOGICAL :: setup_vertical
+
+       setup_vertical = .TRUE.
+       IF (PRESENT(ic_mode))  THEN
+          IF (TRIM(ic_mode) == 'profile')  setup_vertical = .FALSE.
        END IF
 
@@ -1006,5 +985,5 @@
        CASE DEFAULT
 
-          message = "Interpolation mode '" // mode // "' is not supported."
+          message = "Interpolation quantity '" // kind // "' is not supported."
           CALL abort('setup_interpolation', message)
 
@@ -1014,10 +993,8 @@
 
        CALL find_horizontal_neighbours(lat, lon,                               &
-          cosmo_grid % dxi, cosmo_grid % dyi, grid % clat,                     &
-          grid % clon, grid % ii, grid % jj)
+          grid % clat, grid % clon, grid % ii, grid % jj)
 
        CALL compute_horizontal_interp_weights(lat, lon,                        &
-          cosmo_grid % dxi, cosmo_grid % dyi, grid % clat,                     &
-          grid % clon, grid % ii, grid % jj, grid % w_horiz)
+          grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz)
 
 !------------------------------------------------------------------------------
@@ -1044,11 +1021,15 @@
 
        TYPE(grid_definition), POINTER    ::  grid
-       REAL ::  lonmin_pos,lonmax_pos, latmin_pos, latmax_pos
+       REAL(dp)                          ::  lonmin_pos,lonmax_pos, latmin_pos, latmax_pos
+       REAL(dp)                          ::  cosmo_dxi, cosmo_dyi
+
+       cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0))
+       cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0))
 
        ! find horizontal index ranges for profile averaging
-       lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_grid % dxi
-       lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_grid % dxi
-       latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_grid % dyi
-       latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_grid % dyi
+       lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
+       lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
+       latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
+       latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
 
        imin = FLOOR(lonmin_pos)
@@ -1076,33 +1057,10 @@
 ! Description:
 ! ------------
-!> Helper function that computes horizontal domain extend in x or y direction
-!> such that the centres of a boundary grid fall at -dx/2 or lx + dx/2.
-!> 
-!> Input parameters:
-!> -----------------
-!> dxy : grid spacing in x or y direction
-!> lxy : domain length in dxy direction
-!> 
-!> Output parameters:
-!> ------------------
-!> boundary_extent : Domain minimum xymin (maximum xymax) if dxy < 0 (> 0)
-!------------------------------------------------------------------------------!
-    REAL(dp) FUNCTION boundary_extent(dxy, lxy)
-        REAL(dp), INTENT(IN) ::  dxy, lxy
-
-        boundary_extent = 0.5_dp * lxy + SIGN(lxy + ABS(dxy), dxy)
-
-    END FUNCTION boundary_extent
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
 !> Initializes grid_definition-type variables.
 !> 
 !> Input parameters:
 !> -----------------
-!> mode : Initialization mode, distinguishes between PALM-4U and COSMO-DE grids
-!>    as well as grids covering the boundary surfaces. Valid modes are:
+!> kind : Grid kind, distinguishes between PALM-4U and COSMO-DE grids
+!>    as well as grids covering the boundary surfaces. Valid kinds are:
 !>       - 'palm'
 !>       - 'cosmo-de'
@@ -1114,5 +1072,5 @@
 !>    PALM-4U computational domain (i.e. the outer cell faces). The coordinates
 !>    of the generated grid will be inferred from this information taking into
-!>    account the initialization mode. For example, the coordinates of a
+!>    account the initialization mode ic_mode. For example, the coordinates of a
 !>    boundary grid initialized using mode 'eastwest-scalar' will be located in
 !>    planes one half grid point outwards of xmin and xmax.
@@ -1126,12 +1084,13 @@
 !> grid : Grid variable to be initialized.
 !------------------------------------------------------------------------------!
-    SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax, zmin, zmax,  &
-                                    x0, y0, z0, nx, ny, nz, dx, dy, dz, grid, mode)
+    SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax,              &
+                                    x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode)
         CHARACTER(LEN=*), INTENT(IN)           ::  kind
-        CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  mode
+        CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  ic_mode
         INTEGER, INTENT(IN)                    ::  nx, ny, nz
-        REAL(dp), INTENT(IN)                   ::  xmin, xmax, ymin, ymax, zmin, zmax
+        REAL(dp), INTENT(IN)                   ::  xmin, xmax, ymin, ymax
         REAL(dp), INTENT(IN)                   ::  x0, y0, z0
-        REAL(dp), OPTIONAL, INTENT(IN)         ::  dx, dy, dz
+        REAL(dp), INTENT(IN), TARGET, OPTIONAL ::  z(:)
+        REAL(dp), INTENT(IN), TARGET, OPTIONAL ::  zw(:)
         TYPE(grid_definition), INTENT(INOUT)   ::  grid
 
@@ -1142,5 +1101,4 @@
         grid % lx = xmax - xmin
         grid % ly = ymax - ymin
-        grid % lz = zmax - zmin
 
         grid % x0 = x0
@@ -1151,27 +1109,9 @@
 
         CASE('boundary')
-           IF (.NOT. PRESENT(dx))  THEN
-              message = "dx is not present but needed for 'eastwest-scalar' "//&
-                        "initializaton."
+
+           IF (.NOT.PRESENT(z))  THEN
+              message = "z has not been passed but is required for 'boundary' grids"
               CALL abort('init_grid_definition', message)
            END IF
-           IF (.NOT. PRESENT(dy))  THEN
-              message = "dy is not present but needed for 'eastwest-scalar' "//&
-                        "initializaton."
-              CALL abort('init_grid_definition', message)
-           END IF
-           IF (.NOT. PRESENT(dz))  THEN
-              message = "dz is not present but needed for 'eastwest-scalar' "//&
-                        "initializaton."
-              CALL abort('init_grid_definition', message)
-           END IF
-
-           grid % dx  = dx
-           grid % dy  = dy
-           grid % dz  = dz
-
-           grid % dxi = 1.0_dp / grid % dx
-           grid % dyi = 1.0_dp / grid % dy
-           grid % dzi = 1.0_dp / grid % dz
 
            ALLOCATE( grid % x(0:nx) )
@@ -1181,40 +1121,16 @@
            CALL linspace(ymin, ymax, grid % y)
 
-           ALLOCATE( grid % z(0:nz) )
-           CALL linspace(zmin, zmax, grid % z)
+           grid % z => z
 
            ! Allocate neighbour indices and weights
-           IF (TRIM(mode) .NE. 'profile')  THEN
-              ALLOCATE( grid % kk(0:nx, 0:ny, 0:nz, 2) )
+           IF (TRIM(ic_mode) .NE. 'profile')  THEN
+              ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
               grid % kk(:,:,:,:) = -1
 
-              ALLOCATE( grid % w_verti(0:nx, 0:ny, 0:nz, 2) )
+              ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
               grid % w_verti(:,:,:,:) = 0.0_dp
            END IF
         
         CASE('boundary intermediate')
-           IF (.NOT. PRESENT(dx))  THEN
-              message = "dx is not present but needed for 'eastwest-scalar' "//&
-                        "initializaton."
-              CALL abort('init_grid_definition', message)
-           END IF
-           IF (.NOT. PRESENT(dy))  THEN
-              message = "dy is not present but needed for 'eastwest-scalar' "//&
-                        "initializaton."
-              CALL abort('init_grid_definition', message)
-           END IF
-           IF (.NOT. PRESENT(dz))  THEN
-              message = "dz is not present but needed for 'eastwest-scalar' "//&
-                        "initializaton."
-              CALL abort('init_grid_definition', message)
-           END IF
-
-           grid % dx  = dx
-           grid % dy  = dy
-           grid % dz  = dz
-
-           grid % dxi = 1.0_dp / grid % dx
-           grid % dyi = 1.0_dp / grid % dy
-           grid % dzi = 1.0_dp / grid % dz
 
            ALLOCATE( grid % x(0:nx) )
@@ -1224,11 +1140,9 @@
            CALL linspace(ymin, ymax, grid % y)
 
-           ALLOCATE( grid % z(0:nz) )
-           CALL linspace(zmin, zmax, grid % z)
-
            ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny)  )
+
            CALL rotate_to_cosmo(                                               &
-              phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
-              lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude
+              phir = project( grid % y, y0, EARTH_RADIUS ) ,                   & ! = plate-carree latitude
+              lamr = project( grid % x, x0, EARTH_RADIUS ) ,                   & ! = plate-carree longitude
               phip = phi_cn, lamp = lambda_cn,                                 &
               phi  = grid % clat,                                              &
@@ -1249,57 +1163,55 @@
         ! corresponding latitudes and longitudes of the rotated pole grid.
         CASE('palm')
+
+           IF (.NOT.PRESENT(z))  THEN
+              message = "z has not been passed but is required for 'palm' grids"
+              CALL abort('init_grid_definition', message)
+           END IF
+
+           IF (.NOT.PRESENT(zw))  THEN
+              message = "zw has not been passed but is required for 'palm' grids"
+              CALL abort('init_grid_definition', message)
+           END IF
+
            grid % name(1) = 'x and lon'
            grid % name(2) = 'y and lat'
            grid % name(3) = 'z'
 
-           grid % dx  = grid % lx / (nx + 1)
-           grid % dy  = grid % ly / (ny + 1)
-           grid % dz  = grid % lz / (nz + 1)
-
-           grid % dxi = 1.0_dp / grid % dx
-           grid % dyi = 1.0_dp / grid % dy
-           grid % dzi = 1.0_dp / grid % dz
-
-           ALLOCATE( grid % x(0:nx),   grid % y(0:ny),  grid % z(0:nz)  )
-           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny), grid % zw(1:nz) )
-           CALL linspace(xmin + 0.5_dp*grid % dx, xmax - 0.5_dp*grid % dx, grid % x)
-           CALL linspace(ymin + 0.5_dp*grid % dy, ymax - 0.5_dp*grid % dy, grid % y)
-           CALL linspace(zmin + 0.5_dp*grid % dz, zmax - 0.5_dp*grid % dz, grid % z)
-           CALL linspace(xmin + grid % dx, xmax - grid % dx, grid % xu)
-           CALL linspace(ymin + grid % dy, ymax - grid % dy, grid % yv)
-           CALL linspace(zmin + grid % dz, zmax - grid % dz, grid % zw)
+           !TODO: Remove use of global dx, dy, dz variables. Consider
+           !TODO:   associating global x,y, and z arrays.
+           ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
+           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
+           CALL linspace(xmin + 0.5_dp* dx, xmax - 0.5_dp* dx, grid % x)
+           CALL linspace(ymin + 0.5_dp* dy, ymax - 0.5_dp* dy, grid % y)
+           grid % z => z
+           CALL linspace(xmin +  dx, xmax -  dx, grid % xu)
+           CALL linspace(ymin +  dy, ymax -  dy, grid % yv)
+           grid % zw => zw
 
            grid % depths => depths
 
            ! Allocate neighbour indices and weights
-           IF (TRIM(mode) .NE. 'profile')  THEN
-              ALLOCATE( grid % kk(0:nx, 0:ny, 0:nz, 2) )
+           IF (TRIM(ic_mode) .NE. 'profile')  THEN
+              ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
               grid % kk(:,:,:,:) = -1
 
-              ALLOCATE( grid % w_verti(0:nx, 0:ny, 0:nz, 2) )
+              ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
               grid % w_verti(:,:,:,:) = 0.0_dp
            END IF
 
         CASE('palm intermediate')
+
            grid % name(1) = 'x and lon'
            grid % name(2) = 'y and lat'
-           grid % name(3) = 'z'
-
-           grid % dx  = grid % lx / (nx + 1)
-           grid % dy  = grid % ly / (ny + 1)
-           grid % dz  = grid % lz / (nz + 1)
-
-           grid % dxi = 1.0_dp / grid % dx
-           grid % dyi = 1.0_dp / grid % dy
-           grid % dzi = 1.0_dp / grid % dz
-
-           ALLOCATE( grid % x(0:nx),   grid % y(0:ny),  grid % z(0:nz)  )
-           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny), grid % zw(1:nz) )
-           CALL linspace(xmin + 0.5_dp*grid % dx, xmax - 0.5_dp*grid % dx, grid % x)
-           CALL linspace(ymin + 0.5_dp*grid % dy, ymax - 0.5_dp*grid % dy, grid % y)
-           CALL linspace(zmin + 0.5_dp*grid % dz, zmax - 0.5_dp*grid % dz, grid % z)
-           CALL linspace(xmin + grid % dx, xmax - grid % dx, grid % xu)
-           CALL linspace(ymin + grid % dy, ymax - grid % dy, grid % yv)
-           CALL linspace(zmin + grid % dz, zmax - grid % dz, grid % zw)
+           grid % name(3) = 'interpolated hhl or hfl'
+
+           !TODO: Remove use of global dx, dy, dz variables. Consider
+           !TODO:   associating global x,y, and z arrays.
+           ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
+           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
+           CALL linspace(xmin + 0.5_dp*dx, xmax - 0.5_dp*dx, grid % x)
+           CALL linspace(ymin + 0.5_dp*dy, ymax - 0.5_dp*dy, grid % y)
+           CALL linspace(xmin + dx, xmax - dx, grid % xu)
+           CALL linspace(ymin + dy, ymax - dy, grid % yv)
 
            grid % depths => depths
@@ -1355,19 +1267,11 @@
            grid % name(3) = 'height'
 
-           grid % dx  = grid % lx / nx     ! = 0.025 deg, stored in radians
-           grid % dy  = grid % ly / ny     ! = 0.025 deg, stored in radians
-           grid % dz  = 0.0_dp             ! not defined yet
-
-           grid % dxi = 1.0_dp / grid % dx ! [rad^-1]
-           grid % dyi = 1.0_dp / grid % dy ! [rad^-1]
-           grid % dzi = 0.0_dp             ! not defined yet
-
-           ALLOCATE( grid % lon(0:nx),   grid % lat(0:ny),  grid % z(0:nz)  )
-           ALLOCATE( grid % lonu(0:nx),  grid % latv(0:ny), grid % zw(0:nz) )
+           ALLOCATE( grid % lon(0:nx),   grid % lat(0:ny)  )
+           ALLOCATE( grid % lonu(0:nx),  grid % latv(0:ny) )
 
            CALL linspace(xmin, xmax, grid % lon)
            CALL linspace(ymin, ymax, grid % lat)
-           grid % lonu(:) = grid % lon + 0.5_dp * grid % dx
-           grid % latv(:) = grid % lat + 0.5_dp * grid % dy
+           grid % lonu(:) = grid % lon + 0.5_dp * (grid % lx / grid % nx)
+           grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny)
 
            ! Point to heights of half levels (hhl) and compute heights of full
@@ -1384,4 +1288,406 @@
 
     END SUBROUTINE init_grid_definition
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> PALM's stretched vertical grid generator. Forked from PALM revision 3139, see
+!> https://palm.muk.uni-hannover.de/trac/browser/palm/trunk/SOURCE/init_grid.f90?rev=3139
+!> 
+!> This routine computes the levels of scalar points. The levels of the velocity
+!> points are then obtained as the midpoints inbetween using the INIFOR routine
+!> 'modpoints'.
+!------------------------------------------------------------------------------!
+    SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, &
+                           dz_stretch_level_start, dz_stretch_level_end,       &
+                           dz_stretch_factor_array)
+
+       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  z, dz, dz_stretch_factor_array
+       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_level_start, dz_stretch_level_end
+       REAL(dp), INTENT(IN) ::  dz_max, dz_stretch_factor, dz_stretch_level
+
+       INTEGER ::  number_stretch_level_start        !< number of user-specified start levels for stretching
+       INTEGER ::  number_stretch_level_end          !< number of user-specified end levels for stretching
+
+       REAL(dp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end
+       REAL(dp) ::  dz_level_end, dz_stretched
+
+       INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
+       INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
+       INTEGER ::  dz_stretch_level_index = 0
+       INTEGER ::  k, n, number_dz
+!
+!-- Compute height of u-levels from constant grid length and dz stretch factors
+       IF ( dz(1) == -1.0_dp )  THEN
+          message = 'missing dz'
+          CALL abort( 'stretched_z', message)
+       ELSEIF ( dz(1) <= 0.0_dp )  THEN
+          WRITE( message, * ) 'dz=', dz(1),' <= 0.0'
+          CALL abort( 'stretched_z', message)
+       ENDIF
+
+!
+!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
+!-- if it was set by the user
+       IF ( dz_stretch_level /= -9999999.9_dp ) THEN
+          dz_stretch_level_start(1) = dz_stretch_level
+       ENDIF
+       
+!
+!-- Determine number of dz values and stretching levels specified by the
+!-- user to allow right controlling of the stretching mechanism and to
+!-- perform error checks. The additional requirement that dz /= dz_max
+!-- for counting number of user-specified dz values is necessary. Otherwise
+!-- restarts would abort if the old stretching mechanism with dz_stretch_level
+!-- is used (Attention: The user is not allowed to specify a dz value equal
+!-- to the default of dz_max = 999.0). 
+       number_dz = COUNT( dz /= -1.0_dp .AND. dz /= dz_max )
+       number_stretch_level_start = COUNT( dz_stretch_level_start /=           &
+                                           -9999999.9_dp )
+       number_stretch_level_end = COUNT( dz_stretch_level_end /=               &
+                                         9999999.9_dp )
+
+!
+!-- The number of specified end levels +1 has to be the same than the number 
+!-- of specified dz values
+       IF ( number_dz /= number_stretch_level_end + 1 ) THEN
+          WRITE( message, * ) 'The number of values for dz = ',                &
+                              number_dz, 'has to be the same than ',           &
+                              'the number of values for ',                     &
+                              'dz_stretch_level_end + 1 = ',                   &
+                              number_stretch_level_end+1
+          CALL abort( 'stretched_z', message)
+       ENDIF
+    
+!
+!--    The number of specified start levels has to be the same or one less than 
+!--    the number of specified dz values
+       IF ( number_dz /= number_stretch_level_start + 1 .AND.                  &
+            number_dz /= number_stretch_level_start ) THEN
+          WRITE( message, * ) 'The number of values for dz = ',         &
+                              number_dz, 'has to be the same or one ', &
+                              'more than& the number of values for ',  &
+                              'dz_stretch_level_start = ',             &
+                              number_stretch_level_start
+          CALL abort( 'stretched_z', message)
+       ENDIF
+    
+!--    The number of specified start levels has to be the same or one more than 
+!--    the number of specified end levels
+       IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.   &
+            number_stretch_level_start /= number_stretch_level_end ) THEN
+          WRITE( message, * ) 'The number of values for ',              &
+                              'dz_stretch_level_start = ',              &
+                              dz_stretch_level_start, 'has to be the ',&
+                              'same or one more than& the number of ', &
+                              'values for dz_stretch_level_end = ',    &
+                              number_stretch_level_end
+          CALL abort( 'stretched_z', message)
+       ENDIF
+
+!
+!-- Initialize dz for the free atmosphere with the value of dz_max
+       IF ( dz(number_stretch_level_start+1) == -1.0_dp .AND.                     &
+            number_stretch_level_start /= 0 ) THEN 
+          dz(number_stretch_level_start+1) = dz_max
+       ENDIF
+       
+!
+!-- Initialize the stretching factor if (infinitely) stretching in the free 
+!-- atmosphere is desired (dz_stretch_level_end was not specified for the 
+!-- free atmosphere)
+       IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN 
+          dz_stretch_factor_array(number_stretch_level_start) =                   &
+          dz_stretch_factor
+       ENDIF
+
+!-- Allocation of arrays for stretching
+       ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
+
+!
+!--    The stretching region has to be large enough to allow for a smooth
+!--    transition between two different grid spacings
+       DO n = 1, number_stretch_level_start
+          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
+                                        4 * MAX( dz(n),dz(n+1) )
+       ENDDO
+
+       IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
+                 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
+       !IF ( ANY( min_dz_stretch_level_end >      &
+       !          dz_stretch_level_end ) ) THEN
+             message = 'Each dz_stretch_level_end has to be larger '  // &
+                       'than its corresponding value for ' //            &
+                       'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
+                       'to allow for smooth grid stretching'
+             CALL abort('stretched_z', message)
+       ENDIF
+       
+!
+!--    Stretching must not be applied within the prandtl_layer 
+!--    (first two grid points). For the default case dz_stretch_level_start 
+!--    is negative. Therefore the absolut value is checked here.
+       IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_dp ) ) THEN
+          WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',&
+                              'larger than ', dz(1) * 1.5
+             CALL abort( 'stretched_z', message)
+       ENDIF
+
+!
+!--    The stretching has to start and end on a grid level. Therefore 
+!--    user-specified values have to ''interpolate'' to the next lowest level
+       IF ( number_stretch_level_start /= 0 ) THEN
+          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
+                                            dz(1)/2.0) / dz(1) )               &
+                                      * dz(1) + dz(1)/2.0
+       ENDIF
+       
+       IF ( number_stretch_level_start > 1 ) THEN
+          DO n = 2, number_stretch_level_start
+             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
+                                              dz(n) ) * dz(n)
+          ENDDO
+       ENDIF
+       
+       IF ( number_stretch_level_end /= 0 ) THEN
+          DO n = 1, number_stretch_level_end
+             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
+                                            dz(n+1) ) * dz(n+1)
+          ENDDO
+       ENDIF
+ 
+!
+!--    Determine stretching factor if necessary
+       IF ( number_stretch_level_end >= 1 ) THEN 
+          CALL calculate_stretching_factor( number_stretch_level_end, dz,      & 
+                                            dz_stretch_factor,                 &
+                                            dz_stretch_factor_array,           &   
+                                            dz_stretch_level_end,              &
+                                            dz_stretch_level_start )
+       ENDIF
+
+       z(1) = dz(1) * 0.5_dp
+!
+       dz_stretch_level_index = n
+       dz_stretched = dz(1)
+       DO  k = 2, n
+
+          IF ( dz_stretch_level <= z(k-1)  .AND.  dz_stretched < dz_max )  THEN
+
+             dz_stretched = dz_stretched * dz_stretch_factor
+             dz_stretched = MIN( dz_stretched, dz_max )
+
+             IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1
+
+          ENDIF
+
+          z(k) = z(k-1) + dz_stretched
+
+       ENDDO
+!--    Determine u and v height levels considering the possibility of grid
+!--    stretching in several heights.
+       n = 1
+       dz_stretch_level_start_index(:) = UBOUND(z, 1)
+       dz_stretch_level_end_index(:) = UBOUND(z, 1)
+       dz_stretched = dz(1)
+
+!--    The default value of dz_stretch_level_start is negative, thus the first
+!--    condition is always true. Hence, the second condition is necessary.
+       DO  k = 2, UBOUND(z, 1)
+          IF ( dz_stretch_level_start(n) <= z(k-1) .AND.                      &
+               dz_stretch_level_start(n) /= -9999999.9_dp ) THEN
+             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
+             
+             IF ( dz(n) > dz(n+1) ) THEN
+                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
+             ELSE
+                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
+             ENDIF
+             
+             IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) )            &
+             dz_stretch_level_start_index(n) = k-1
+             
+          ENDIF
+          
+          z(k) = z(k-1) + dz_stretched
+          
+!
+!--       Make sure that the stretching ends exactly at dz_stretch_level_end 
+          dz_level_end = ABS( z(k) - dz_stretch_level_end(n) ) 
+          
+          IF ( dz_level_end < dz(n+1)/3.0 ) THEN
+             z(k) = dz_stretch_level_end(n)
+             dz_stretched = dz(n+1)
+             dz_stretch_level_end_index(n) = k
+             n = n + 1             
+          ENDIF
+       ENDDO
+
+       DEALLOCATE( min_dz_stretch_level_end )
+
+    END SUBROUTINE stretched_z
+
+
+! Description: [PALM subroutine]
+! -----------------------------------------------------------------------------!
+!> Calculation of the stretching factor through an iterative method. Ideas were 
+!> taken from the paper "Regional stretched grid generation and its application
+!> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
+!> system of equations has two variables (r,l) but four requirements 
+!> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
+!> results into an overdetermined system. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE calculate_stretching_factor( number_end, dz, dz_stretch_factor,    &
+                                         dz_stretch_factor_array,              &   
+                                         dz_stretch_level_end,                 &
+                                         dz_stretch_level_start )
+ 
+    REAL(dp), DIMENSION(:), INTENT(IN)    ::  dz
+    REAL(dp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_factor_array
+    REAL(dp), DIMENSION(:), INTENT(IN)    ::  dz_stretch_level_end, dz_stretch_level_start
+    REAL(dp)                              ::  dz_stretch_factor
+  
+    INTEGER ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached  
+    INTEGER ::  l_rounded   !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2 
+    INTEGER ::  n           !< loop variable for stretching
+    
+    INTEGER, INTENT(IN) ::  number_end !< number of user-specified end levels for stretching
+        
+    REAL(dp) ::  delta_l               !< absolute difference between l and l_rounded
+    REAL(dp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
+    REAL(dp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible) 
+    REAL(dp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration 
+    REAL(dp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
+    REAL(dp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
+    REAL(dp) ::  numerator             !< numerator of the quotient
+    REAL(dp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
+    REAL(dp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
+    
+    REAL(dp) ::  dz_stretch_factor_array_2(9) = 1.08_dp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1 
+    
+    REAL(dp), PARAMETER ::  stretch_factor_interval = 1.0E-06  !< interval for sampling possible stretching factors
+    REAL(dp), PARAMETER ::  stretch_factor_lower_limit = 0.88  !< lowest possible stretching factor
+    REAL(dp), PARAMETER ::  stretch_factor_upper_limit = 1.12  !< highest possible stretching factor
+ 
+ 
+    l = 0
+    DO  n = 1, number_end
+    
+       iterations = 1
+       stretch_factor_1 = 1.0 
+       stretch_factor_2 = 1.0
+       delta_total_old = 1.0
+       
+       IF ( dz(n) > dz(n+1) ) THEN
+          DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 
+             
+             stretch_factor_1 = 1.0 - iterations * stretch_factor_interval
+             distance = ABS( dz_stretch_level_end(n) -                   &
+                        dz_stretch_level_start(n) ) 
+             numerator = distance*stretch_factor_1/dz(n) +               &
+                         stretch_factor_1 - distance/dz(n)
+             
+             IF ( numerator > 0.0 ) THEN
+                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
+                l_rounded = NINT( l )
+                delta_l = ABS( l_rounded - l ) / l
+             ENDIF
+             
+             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
+             
+             delta_stretch_factor = ABS( stretch_factor_1 -              &
+                                         stretch_factor_2 ) /            &
+                                    stretch_factor_2
+             
+             delta_total_new = delta_l + delta_stretch_factor
+
+!
+!--                stretch_factor_1 is taken to guarantee that the stretching
+!--                procedure ends as close as possible to dz_stretch_level_end.
+!--                stretch_factor_2 would guarantee that the stretched dz(n) is
+!--                equal to dz(n+1) after l_rounded grid levels.
+             IF (delta_total_new < delta_total_old) THEN
+                dz_stretch_factor_array(n) = stretch_factor_1
+                dz_stretch_factor_array_2(n) = stretch_factor_2
+                delta_total_old = delta_total_new
+             ENDIF
+             
+             iterations = iterations + 1
+            
+          ENDDO
+             
+       ELSEIF ( dz(n) < dz(n+1) ) THEN 
+          DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
+                     
+             stretch_factor_1 = 1.0 + iterations * stretch_factor_interval
+             distance = ABS( dz_stretch_level_end(n) -                      &
+                        dz_stretch_level_start(n) ) 
+             numerator = distance*stretch_factor_1/dz(n) +                  &
+                         stretch_factor_1 - distance/dz(n)
+             
+             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
+             l_rounded = NINT( l )
+             delta_l = ABS( l_rounded - l ) / l
+             
+             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
+
+             delta_stretch_factor = ABS( stretch_factor_1 -                 &
+                                        stretch_factor_2 ) /                &
+                                        stretch_factor_2
+             
+             delta_total_new = delta_l + delta_stretch_factor
+             
+!
+!--                stretch_factor_1 is taken to guarantee that the stretching
+!--                procedure ends as close as possible to dz_stretch_level_end.
+!--                stretch_factor_2 would guarantee that the stretched dz(n) is
+!--                equal to dz(n+1) after l_rounded grid levels.
+             IF (delta_total_new < delta_total_old) THEN
+                dz_stretch_factor_array(n) = stretch_factor_1
+                dz_stretch_factor_array_2(n) = stretch_factor_2
+                delta_total_old = delta_total_new
+             ENDIF
+             
+             iterations = iterations + 1
+          ENDDO
+          
+       ELSE
+          message = 'Two adjacent values of dz must be different'
+          CALL abort( 'calculate_stretching_factor', message)
+       ENDIF
+
+!
+!--    Check if also the second stretching factor fits into the allowed
+!--    interval. If not, print a warning for the user.
+       IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR.     & 
+            dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN
+          WRITE( message, * ) 'stretch_factor_2 = ',                    &
+                                     dz_stretch_factor_array_2(n), ' which is',&
+                                     ' responsible for exactly reaching& dz =',&
+                                      dz(n+1), 'after a specific amount of',   & 
+                                     ' grid levels& exceeds the upper',        &
+                                     ' limit =', stretch_factor_upper_limit,   &
+                                     ' &or lower limit = ',                    &
+                                     stretch_factor_lower_limit
+          CALL abort( 'calculate_stretching_factor', message )
+            
+       ENDIF
+    ENDDO
+        
+ END SUBROUTINE calculate_stretching_factor
+
+    SUBROUTINE midpoints(z, zw)
+
+        REAL(dp), INTENT(IN)  ::  z(0:)
+        REAL(dp), INTENT(OUT) ::  zw(1:)
+
+        INTEGER ::  k
+
+        DO k = 1, UBOUND(zw, 1)
+           zw(k) = 0.5_dp * (z(k-1) + z(k))
+        END DO
+
+    END SUBROUTINE midpoints
 
 
@@ -1409,16 +1715,17 @@
        )
 
-       !potential temperature, surface pressure
+       !potential temperature, surface pressure, including nudging and subsidence
        io_group_list(3) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          out_vars = [output_var_table(3:8), output_var_table(42:42)],         &
+          out_vars = [output_var_table(3:8), output_var_table(42:42),          &
+                      output_var_table(49:51)],                                &
           in_var_list = (/input_var_table(3), input_var_table(17)/),           &
           kind = 'temperature'                                                 &
        )
 
-       !specific humidity
+       !specific humidity including nudging and subsidence
        io_group_list(4) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          out_vars = output_var_table(9:14),                                   &
+          out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
           in_var_list = input_var_table(4:4),                                  &
           kind = 'scalar'                                                      &
@@ -1428,5 +1735,5 @@
        io_group_list(5) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          out_vars = [output_var_table(15:26), output_var_table(43:44)],       &
+          out_vars = [output_var_table(15:26), output_var_table(43:46)],       &
           !out_vars = output_var_table(15:20),                                  &
           in_var_list = input_var_table(5:6),                                  &
@@ -1444,8 +1751,8 @@
        !io_group_list(6) % to_be_processed = .FALSE.
    
-       !w velocity
+       !w velocity and subsidence and w nudging
        io_group_list(7) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          out_vars = output_var_table(27:32),                                  &
+          out_vars = [output_var_table(27:32), output_var_table(47:48)],       &
           in_var_list = input_var_table(7:7),                                  &
           kind = 'scalar'                                                      &
@@ -1459,4 +1766,5 @@
           kind = 'accumulated'                                                 &
        )
+       io_group_list(8) % to_be_processed = .FALSE.
 
        !snow
@@ -1467,4 +1775,5 @@
           kind = 'accumulated'                                                 &
        )
+       io_group_list(9) % to_be_processed = .FALSE.
 
        !graupel
@@ -1475,4 +1784,5 @@
           kind = 'accumulated'                                                 &
        )
+       io_group_list(10) % to_be_processed = .FALSE.
 
        !evapotranspiration
@@ -1483,4 +1793,5 @@
           kind = 'accumulated'                                                 &
        )
+       io_group_list(11) % to_be_processed = .FALSE.
 
        !2m air temperature
@@ -1517,4 +1828,5 @@
           kind = 'running average'                                             &
        )
+       io_group_list(15) % to_be_processed = .FALSE.
 
        !lw radiation balance
@@ -1525,4 +1837,5 @@
           kind = 'running average'                                             &
        )
+       io_group_list(16) % to_be_processed = .FALSE.
 
     END SUBROUTINE setup_io_groups
@@ -1562,4 +1875,6 @@
        CALL report('fini_grids', 'Deallocating grids')
        
+       DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
+
        DEALLOCATE(palm_grid%x,  palm_grid%y,  palm_grid%z,                     &
                   palm_grid%xu, palm_grid%yv, palm_grid%zw,                    &
@@ -1572,6 +1887,6 @@
                   palm_intermediate%clonu, palm_intermediate%clatu)
 
-       DEALLOCATE(cosmo_grid%lon,  cosmo_grid%lat,  cosmo_grid%z,              &
-                  cosmo_grid%lonu, cosmo_grid%latv, cosmo_grid%zw,             &
+       DEALLOCATE(cosmo_grid%lon,  cosmo_grid%lat,                             &
+                  cosmo_grid%lonu, cosmo_grid%latv,                            &
                   cosmo_grid%hfl)
 
@@ -1584,17 +1899,17 @@
 !> Initializes the the variable list.
 !------------------------------------------------------------------------------!
-    SUBROUTINE setup_variable_tables(mode)
-       CHARACTER(LEN=*), INTENT(IN) ::  mode
+    SUBROUTINE setup_variable_tables(ic_mode)
+       CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
        TYPE(nc_var), POINTER        ::  var
 
-       IF (TRIM(start_date) == '')  THEN
+       IF (TRIM(cfg % start_date) == '')  THEN
           message = 'Simulation start date has not been set.'
           CALL abort('setup_variable_tables', message)
        END IF
 
-       nc_source_text = 'COSMO-DE analysis from ' // TRIM(start_date)
+       nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg % start_date)
 
        n_invar = 17
-       n_outvar = 44
+       n_outvar = 55
        ALLOCATE( input_var_table(n_invar) )
        ALLOCATE( output_var_table(n_outvar) )
@@ -1693,4 +2008,8 @@
 !- Section 2: NetCDF output variables
 !------------------------------------------------------------------------------
+!
+!------------------------------------------------------------------------------
+! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions
+!------------------------------------------------------------------------------
        output_var_table(1) = init_nc_var(                                      &
           name              = 'init_soil_t',                                   &
@@ -1718,5 +2037,5 @@
 
        output_var_table(3) = init_nc_var(                                      &
-          name              = 'init_pt',                                       &
+          name              = 'init_atmosphere_pt',                            &
           std_name          = "",                                              &
           long_name         = "initial potential temperature",                 &
@@ -1727,7 +2046,7 @@
           grid              = palm_grid,                                       &
           intermediate_grid = palm_intermediate,                               &
-          is_profile = (TRIM(mode) == 'profile')                               &
-       )
-       IF (TRIM(mode) == 'profile')  THEN
+          is_profile = (TRIM(ic_mode) == 'profile')                            &
+       )
+       IF (TRIM(ic_mode) == 'profile')  THEN
           output_var_table(3) % grid => scalar_profile_grid
           output_var_table(3) % intermediate_grid => scalar_profile_intermediate
@@ -1795,5 +2114,5 @@
 
        output_var_table(9) = init_nc_var(                                      &
-          name              = 'init_qv',                                       &
+          name              = 'init_atmosphere_qv',                            &
           std_name          = "",                                              &
           long_name         = "initial specific humidity",                     &
@@ -1804,7 +2123,7 @@
           grid              = palm_grid,                                       &
           intermediate_grid = palm_intermediate,                               &
-          is_profile = (TRIM(mode) == 'profile')                               &
-       )
-       IF (TRIM(mode) == 'profile')  THEN
+          is_profile = (TRIM(ic_mode) == 'profile')                            &
+       )
+       IF (TRIM(ic_mode) == 'profile')  THEN
           output_var_table(9) % grid => scalar_profile_grid
           output_var_table(9) % intermediate_grid => scalar_profile_intermediate
@@ -1872,5 +2191,5 @@
 
        output_var_table(15) = init_nc_var(                                     &
-          name              = 'init_u',                                        &
+          name              = 'init_atmosphere_u',                             &
           std_name          = "",                                              &
           long_name         = "initial wind component in x direction",         &
@@ -1881,7 +2200,7 @@
           grid              = u_initial_grid,                                  &
           intermediate_grid = u_initial_intermediate,                          &
-          is_profile = (TRIM(mode) == 'profile')                               &
-       )
-       IF (TRIM(mode) == 'profile')  THEN
+          is_profile = (TRIM(ic_mode) == 'profile')                            &
+       )
+       IF (TRIM(ic_mode) == 'profile')  THEN
           output_var_table(15) % grid => scalar_profile_grid
           output_var_table(15) % intermediate_grid => scalar_profile_intermediate
@@ -1949,5 +2268,5 @@
 
        output_var_table(21) = init_nc_var(                                     &
-          name              = 'init_v',                                        &
+          name              = 'init_atmosphere_v',                             &
           std_name          = "",                                              &
           long_name         = "initial wind component in y direction",         &
@@ -1958,7 +2277,7 @@
           grid              = v_initial_grid,                                  &
           intermediate_grid = v_initial_intermediate,                          &
-          is_profile = (TRIM(mode) == 'profile')                               &
-       )
-       IF (TRIM(mode) == 'profile')  THEN
+          is_profile = (TRIM(ic_mode) == 'profile')                            &
+       )
+       IF (TRIM(ic_mode) == 'profile')  THEN
           output_var_table(21) % grid => scalar_profile_grid
           output_var_table(21) % intermediate_grid => scalar_profile_intermediate
@@ -2026,5 +2345,5 @@
 
        output_var_table(27) = init_nc_var(                                     &
-          name              = 'init_w',                                        &
+          name              = 'init_atmosphere_w',                             &
           std_name          = "",                                              &
           long_name         = "initial wind component in z direction",         &
@@ -2035,7 +2354,7 @@
           grid              = w_initial_grid,                                  &
           intermediate_grid = w_initial_intermediate,                          &
-          is_profile = (TRIM(mode) == 'profile')                               &
-       )
-       IF (TRIM(mode) == 'profile')  THEN
+          is_profile = (TRIM(ic_mode) == 'profile')                            &
+       )
+       IF (TRIM(ic_mode) == 'profile')  THEN
           output_var_table(27) % grid => w_profile_grid
           output_var_table(27) % intermediate_grid => w_profile_intermediate
@@ -2209,5 +2528,8 @@
           intermediate_grid = palm_intermediate                                &
        )
-
+!
+!------------------------------------------------------------------------------
+! Section 2.2: Idealized large-scale forcings
+!------------------------------------------------------------------------------
        output_var_table(42) = init_nc_var(                                     &
           name              = 'surface_forcing_surface_pressure',              &
@@ -2227,9 +2549,9 @@
           long_name         = "geostrophic wind (u component)",                &
           units             = "m/s",                                           &
-          kind              = "profile",                                       &
-          input_id          = 1,                                               &
-          output_file       = output_file,                                     &
-          grid              = palm_grid,                                       &
-          intermediate_grid = palm_intermediate                                &
+          kind              = "constant scalar profile",                       &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
        )
 
@@ -2239,13 +2561,160 @@
           long_name         = "geostrophic wind (v component)",                &
           units             = "m/s",                                           &
-          kind              = "profile",                                       &
-          input_id          = 1,                                               &
-          output_file       = output_file,                                     &
-          grid              = palm_grid,                                       &
-          intermediate_grid = palm_intermediate                                &
-       )
+          kind              = "constant scalar profile",                       &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+
+       output_var_table(45) = init_nc_var(                                     &
+          name              = 'nudging_u',                                     &
+          std_name          = "",                                              &
+          long_name         = "wind component in x direction",                 &
+          units             = "m/s",                                           &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(45) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(46) = init_nc_var(                                     &
+          name              = 'nudging_v',                                     &
+          std_name          = "",                                              &
+          long_name         = "wind component in y direction",                 &
+          units             = "m/s",                                           &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(46) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(47) = init_nc_var(                                     &
+          name              = 'ls_forcing_sub_w',                              &
+          std_name          = "",                                              &
+          long_name         = "subsidence velocity of w",                      &
+          units             = "m/s",                                           &
+          kind              = "large-scale w forcing",                         &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = w_profile_grid,                                  &
+          intermediate_grid = w_profile_intermediate                           &
+       )
+       output_var_table(47) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(48) = init_nc_var(                                     &
+          name              = 'nudging_w',                                     &
+          std_name          = "",                                              &
+          long_name         = "wind component in w direction",                 &
+          units             = "m/s",                                           &
+          kind              = "large-scale w forcing",                         &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = w_profile_grid,                                  &
+          intermediate_grid = w_profile_intermediate                           &
+       )
+       output_var_table(48) % to_be_processed = ls_forcing_variables_required
+
+
+       output_var_table(49) = init_nc_var(                                     &
+          name              = 'ls_forcing_adv_pt',                             &
+          std_name          = "",                                              &
+          long_name         = "advection of potential temperature",            &
+          units             = "K/s",                                           &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(49) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(50) = init_nc_var(                                     &
+          name              = 'ls_forcing_sub_pt',                             &
+          std_name          = "",                                              &
+          long_name         = "subsidence velocity of potential temperature",  &
+          units             = "K/s",                                           &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(50) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(51) = init_nc_var(                                     &
+          name              = 'nudging_pt',                                    &
+          std_name          = "",                                              &
+          long_name         = "potential temperature",                         &
+          units             = "K",                                             &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(51) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(52) = init_nc_var(                                     &
+          name              = 'ls_forcing_adv_qv',                             &
+          std_name          = "",                                              &
+          long_name         = "advection of specific humidity",                &
+          units             = "kg/kg/s",                                       &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(52) % to_be_processed = ls_forcing_variables_required
+
+
+       output_var_table(53) = init_nc_var(                                     &
+          name              = 'ls_forcing_sub_qv',                             &
+          std_name          = "",                                              &
+          long_name         = "subsidence velocity of specific humidity",      &
+          units             = "kg/kg/s",                                       &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(53) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(54) = init_nc_var(                                     &
+          name              = 'nudging_qv',                                    &
+          std_name          = "",                                              &
+          long_name         = "specific humidity",                             &
+          units             = "kg/kg",                                         &
+          kind              = "large-scale scalar forcing",                    &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(54) % to_be_processed = ls_forcing_variables_required
+
+       output_var_table(55) = init_nc_var(                                     &
+          name              = 'nudging_tau',                                   &
+          std_name          = "",                                              &
+          long_name         = "nudging relaxation time scale",                 &
+          units             = "s",                                             &
+          kind              = "constant scalar profile",                       &
+          input_id          = 1,                                               &
+          output_file       = output_file,                                     &
+          grid              = scalar_profile_grid,                             &
+          intermediate_grid = scalar_profile_intermediate                      &
+       )
+       output_var_table(55) % to_be_processed = ls_forcing_variables_required
+
 
        ! Attributes shared among all variables
        output_var_table(:) % source = nc_source_text
+
 
     END SUBROUTINE setup_variable_tables
@@ -2260,5 +2729,6 @@
 !------------------------------------------------------------------------------!
     FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id,     &
-                         grid, intermediate_grid, output_file, is_profile) RESULT(var)
+                         grid, intermediate_grid, output_file, is_profile      &
+             ) RESULT(var)
 
        CHARACTER(LEN=*), INTENT(IN)      ::  name, std_name, long_name, units, kind
@@ -2367,6 +2837,6 @@
           var % task            = "average profile"
 
-       CASE ( 'surface forcing' )
-          var % lod             = 1
+       CASE( 'surface forcing' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time
@@ -2377,6 +2847,6 @@
           var % task            = "interpolate_2d"
 
-       CASE ( 'left scalar', 'right scalar') ! same as right
-          var % lod             = 2
+       CASE( 'left scalar', 'right scalar') ! same as right
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time
@@ -2389,6 +2859,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'north scalar', 'south scalar') ! same as south
-          var % lod             = 2
+       CASE( 'north scalar', 'south scalar') ! same as south
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time
@@ -2401,6 +2871,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'top scalar', 'top w' )
-          var % lod             = 2
+       CASE( 'top scalar', 'top w' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time
@@ -2413,6 +2883,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'left u', 'right u' )
-          var % lod             = 2
+       CASE( 'left u', 'right u' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time
@@ -2425,6 +2895,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'north u', 'south u' )
-          var % lod             = 2
+       CASE( 'north u', 'south u' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time    !t
@@ -2437,6 +2907,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'top u' )
-          var % lod             = 2
+       CASE( 'top u' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time    !t
@@ -2449,6 +2919,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'left v', 'right v' )
-          var % lod             = 2
+       CASE( 'left v', 'right v' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time
@@ -2461,6 +2931,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'north v', 'south v' )
-          var % lod             = 2
+       CASE( 'north v', 'south v' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time    !t
@@ -2473,6 +2943,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'top v' )
-          var % lod             = 2
+       CASE( 'top v' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time    !t
@@ -2485,6 +2955,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'left w', 'right w')
-          var % lod             = 2
+       CASE( 'left w', 'right w')
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time
@@ -2497,6 +2967,6 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'north w', 'south w' )
-          var % lod             = 2
+       CASE( 'north w', 'south w' )
+          var % lod             = -1
           var % ndim            = 3
           var % dimids(3)       = output_file % dimid_time    !t
@@ -2509,5 +2979,5 @@
           var % task            = "interpolate_3d"
 
-       CASE ( 'time series' )
+       CASE( 'time series' )
           var % lod             = 0
           var % ndim            = 1
@@ -2517,6 +2987,6 @@
           var % task            = "average scalar"
 
-       CASE ( 'profile' )
-          var % lod             = 2
+       CASE( 'constant scalar profile' )
+          var % lod             = -1
           var % ndim            = 2
           var % dimids(2)       = output_file % dimid_time    !t
@@ -2525,5 +2995,25 @@
           var % dimvarids(1)    = output_file % dimvarids_scl(3)
           var % to_be_processed = .TRUE.
-          var % task            = "profile"
+          var % task            = "set profile"
+
+       CASE( 'large-scale scalar forcing' )
+          var % lod             = -1
+          var % ndim            = 2
+          var % dimids(2)       = output_file % dimid_time    !t
+          var % dimids(1)       = output_file % dimids_scl(3) !z
+          var % dimvarids(2)    = output_file % dimvarid_time
+          var % dimvarids(1)    = output_file % dimvarids_scl(3)
+          var % to_be_processed = ls_forcing_variables_required
+          var % task            = "average large-scale profile"
+
+       CASE( 'large-scale w forcing' )
+          var % lod             = -1
+          var % ndim            = 2
+          var % dimids(2)       = output_file % dimid_time    !t
+          var % dimids(1)       = output_file % dimids_vel(3) !z
+          var % dimvarids(2)    = output_file % dimvarid_time
+          var % dimvarids(1)    = output_file % dimvarids_vel(3)
+          var % to_be_processed = ls_forcing_variables_required
+          var % task            = "average large-scale profile"
 
        CASE DEFAULT
@@ -2560,6 +3050,6 @@
 
 
-    SUBROUTINE input_file_list(start_date_string, start_hour, end_hour,        &
-                               step_hour, path, prefix, suffix, file_list)
+    SUBROUTINE get_input_file_list(start_date_string, start_hour, end_hour,        &
+                                   step_hour, path, prefix, suffix, file_list)
 
        CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
@@ -2568,22 +3058,21 @@
        CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
 
-       INTEGER             ::  number_of_files, hour, i
+       INTEGER             ::  number_of_intervals, hour, i
        CHARACTER(LEN=DATE) ::  date_string
 
-       number_of_files = end_hour - start_hour + 1
-
-       ALLOCATE( file_list(number_of_files) )
-
-       DO i = 1, number_of_files
-          hour = start_hour + (i-1) * step_hour
+       number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
+       ALLOCATE( file_list(number_of_intervals + 1) )
+
+       DO i = 0, number_of_intervals
+          hour = start_hour + i * step_hour
           date_string = add_hours_to(start_date_string, hour)
 
-          file_list(i) = TRIM(path) // TRIM(prefix) // TRIM(date_string) //    &
-                         TRIM(suffix) // '.nc'
-          message = "Set up input file name '" // TRIM(file_list(i)) //"'"
+          file_list(i+1) = TRIM(path) // TRIM(prefix) // TRIM(date_string) //    &
+                           TRIM(suffix) // '.nc'
+          message = "Set up input file name '" // TRIM(file_list(i+1)) //"'"
           CALL report('input_file_list', message)
        END DO
 
-    END SUBROUTINE input_file_list
+    END SUBROUTINE get_input_file_list
 
 
@@ -2607,5 +3096,5 @@
        
        REAL(dp), ALLOCATABLE                       ::  basic_state_pressure(:)
-       TYPE(container), ALLOCATABLE                ::  compute_buffer(:)
+       TYPE(container), ALLOCATABLE                ::  preprocess_buffer(:)
        INTEGER                                     ::  hour, dt
        INTEGER                                     ::  i, j, k
@@ -2617,6 +3106,6 @@
        
        CASE( 'velocities' )
-          ! Allocate a compute puffer with the same number of arrays as the input
-          ALLOCATE( compute_buffer( SIZE(input_buffer) ) )
+          ! Allocate a compute buffer with the same number of arrays as the input
+          ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) )
 
           ! Allocate u and v arrays with scalar dimensions
@@ -2624,6 +3113,6 @@
           ny = SIZE(input_buffer(1) % array, 2)
           nz = SIZE(input_buffer(1) % array, 3)
-          ALLOCATE( compute_buffer(1) % array(nx, ny, nz) ) ! u buffer
-          ALLOCATE( compute_buffer(2) % array(nx, ny, nz) ) ! v buffer
+          ALLOCATE( preprocess_buffer(1) % array(nx, ny, nz) ) ! u buffer
+          ALLOCATE( preprocess_buffer(2) % array(nx, ny, nz) ) ! v buffer
 
  CALL run_control('time', 'alloc')
@@ -2632,23 +3121,36 @@
           CALL centre_velocities( u_face = input_buffer(1) % array,            &
                                   v_face = input_buffer(2) % array,            &
-                                  u_centre = compute_buffer(1) % array,        &
-                                  v_centre = compute_buffer(2) % array )
+                                  u_centre = preprocess_buffer(1) % array,     &
+                                  v_centre = preprocess_buffer(2) % array )
           
-          ! rotate U and V to PALM-4U orientation and overwrite U and V with
-          ! rotated velocities
-          DO k = 1, nz
-          DO j = 2, ny
-          DO i = 2, nx
-             CALL uv2uvrot( urot = compute_buffer(1) % array(i,j,k),           &
-                            vrot = compute_buffer(2) % array(i,j,k),           &
-                            rlat = cosmo_grid % lat(j-1),                      &
-                            rlon = cosmo_grid % lon(i-1),                      &
-                            pollat = phi_cn,                                   &
-                            pollon = lambda_cn,                                &
-                            u = input_buffer(1) % array(i,j,k),                &
-                            v = input_buffer(2) % array(i,j,k) )
-          END DO
-          END DO
-          END DO
+          cfg % rotation_method = 'rotated-pole'
+          SELECT CASE(cfg % rotation_method)
+
+          CASE('rotated-pole')
+             ! rotate U and V to PALM-4U orientation and overwrite U and V with
+             ! rotated velocities
+             DO k = 1, nz
+             DO j = 2, ny
+             DO i = 2, nx
+                CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k),     &
+                               vrot = preprocess_buffer(2) % array(i,j,k),     &
+                               rlat = cosmo_grid % lat(j-1),                   &
+                               rlon = cosmo_grid % lon(i-1),                   &
+                               pollat = phi_cn,                                &
+                               pollon = lambda_cn,                             &
+                               u = input_buffer(1) % array(i,j,k),             &
+                               v = input_buffer(2) % array(i,j,k) )
+             END DO
+             END DO
+             END DO
+
+          CASE DEFAULT
+             message = "Rotation method '" // TRIM(cfg % rotation_method) //   &
+                "' not recognized."
+             CALL abort('preprocess', message)
+
+          END SELECT
+
+          ! set values 
           input_buffer(1) % array(1,:,:) = 0.0_dp
           input_buffer(2) % array(1,:,:) = 0.0_dp
@@ -2659,5 +3161,5 @@
  CALL run_control('time', 'comp')
 
-          DEALLOCATE( compute_buffer )
+          DEALLOCATE( preprocess_buffer )
  CALL run_control('time', 'alloc')
 
Index: palm/trunk/UTIL/inifor/src/inifor.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/inifor.f90	(revision 3182)
@@ -21,4 +21,7 @@
 ! Current revisions:
 ! -----------------
+! Introduced new PALM grid stretching
+! Renamend initial-condition mode variable 'mode' to 'ic_mode'
+! Improved log messages
 ! 
 ! 
@@ -48,8 +51,8 @@
         ONLY:  setup_parameters, setup_grids, setup_variable_tables,           &
                setup_io_groups, fini_grids, fini_variables, fini_io_groups,    &
-               fini_file_lists, preprocess,                                    &
+               fini_file_lists, preprocess, origin_lon, origin_lat,            &
                output_file, io_group_list, output_var_table,                   &
-               cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, mode,            &
-               imin, imax, jmin, jmax
+               cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, cfg,             &
+               average_imin, average_imax, average_jmin, average_jmax
 
     USE io
@@ -84,10 +87,11 @@
 
     ! Initialize the netCDF output file and define dimensions
-    CALL setup_netcdf_dimensions(output_file, palm_grid)
+    CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date,    &
+                                 origin_lon, origin_lat)
  CALL run_control('time', 'write')
 
     ! Set up the tables containing the input and output variables and set
     ! the corresponding netCDF dimensions for each output variable
-    CALL setup_variable_tables(mode)
+    CALL setup_variable_tables(cfg % ic_mode)
  CALL run_control('time', 'write')
 
@@ -95,5 +99,5 @@
     CALL setup_netcdf_variables(output_file % name, output_var_table)
 
-    CALL setup_io_groups() 
+    CALL setup_io_groups()
  CALL run_control('time', 'init')
 
@@ -118,4 +122,5 @@
  CALL run_control('time', 'comp')
 
+             !TODO: move this assertion into 'preprocess'.
              IF ( .NOT. ALL(input_buffer(:) % is_preprocessed .AND. .TRUE.) )  THEN
                 message = "Input buffers for group '" // TRIM(group % kind) // &
@@ -159,5 +164,6 @@
                       CASE DEFAULT
 
-                          CALL abort("main loop", 'Not a soil variable')
+                          message = "'" // TRIM(output_var % kind) // "' is not a soil variable"
+                          CALL abort("main loop", message)
 
                       END SELECT
@@ -173,5 +179,5 @@
                       ALLOCATE( output_arr( 0:output_var % grid % nx,          &
                                             0:output_var % grid % ny,          &
-                                            0:output_var % grid % nz ) )
+                                            1:output_var % grid % nz ) )
 
  CALL run_control('time', 'alloc')
@@ -187,5 +193,5 @@
                       ALLOCATE( output_arr( 0:output_var % grid % nx,          &
                                             0:output_var % grid % ny,          &
-                                            0:output_var % grid % nz ) )
+                                            1:output_var % grid % nz ) )
  CALL run_control('time', 'alloc')
                       
@@ -193,5 +199,6 @@
                       CALL average_profile(                                    &
                          input_buffer(output_var % input_id) % array(:,:,:),   &
-                         output_arr(:,:,:), imin, imax, jmin, jmax,            &
+                         output_arr(:,:,:), average_imin, average_imax,        &
+                         average_jmin, average_jmax,                           &
                          output_var % intermediate_grid,                       &
                          output_var % grid)
@@ -205,7 +212,7 @@
  CALL run_control('time', 'comp')
 
-                   CASE ( 'profile' )
+                   CASE ( 'set profile' )
                       
-                      ALLOCATE( output_arr( 1, 1, 0:nz ) )
+                      ALLOCATE( output_arr( 1, 1, 1:nz ) )
  CALL run_control('time', 'alloc')
 
@@ -217,4 +224,7 @@
                       CASE('ls_forcing_vg')
                           output_arr(1, 1, :) = vg
+
+                      CASE('nudging_tau')
+                          output_arr(1, 1, :) = NUDGING_TAU
 
                       CASE DEFAULT
@@ -225,4 +235,10 @@
                       END SELECT
  CALL run_control('time', 'comp')
+
+                   CASE('average large-scale profile')
+                      message = "Averaging of large-scale forcing profiles " //&
+                                "has not been implemented, yet."
+                      CALL abort('main loop', message)
+                      !ALLOCATE( output_arr( 1, 1, 1:nz ) )
 
                    CASE DEFAULT
@@ -269,5 +285,5 @@
        ELSE 
 
-          message = "Skipping IO group '" // TRIM(group % kind) // "'"
+          message = "Skipping IO group " // TRIM(str(igroup)) // " '" // TRIM(group % kind) // "'"
           IF ( ALLOCATED(group % in_var_list) )  THEN
               message = TRIM(message) // " with input variable '" //           &
@@ -291,5 +307,5 @@
  CALL run_control('report', 'void')
 
-    message = "Finished writing forcing file '" // TRIM(output_file % name) // &
+    message = "Finished writing dynamic driver '" // TRIM(output_file % name) // &
               "' successfully."
     CALL report('main loop', message)
Index: palm/trunk/UTIL/inifor/src/io.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/io.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/io.f90	(revision 3182)
@@ -21,6 +21,17 @@
 ! Current revisions:
 ! -----------------
+! Introduced new PALM grid stretching
+! Updated variable names and metadata for PIDS v1.9 compatibility
+! Improved handling of the start date string
+! Better compatibility with older Intel compilers:
+! - avoiding implicit array allocation with new get_netcdf_variable()
+!   subroutine instead of function
+! Improved command line interface:
+! - Added configuration validation
+! - New options to configure input file prefixes
+! - GNU-style short and long option names
+! - Added version and copyright output
 ! 
-! 
+!
 ! Former revisions:
 ! -----------------
@@ -43,13 +54,42 @@
     USE control
     USE defs,                                                                  &
-        ONLY:  DATE, SNAME, PATH, PI, dp, TO_RADIANS, TO_DEGREES, VERSION
+        ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
     USE netcdf
     USE types
     USE util,                                                                  &
-        ONLY:  reverse, str
+        ONLY:  reverse, str, real_to_str
 
     IMPLICIT NONE
 
+    INTERFACE get_netcdf_variable
+        MODULE PROCEDURE get_netcdf_variable_int
+        MODULE PROCEDURE get_netcdf_variable_real
+    END INTERFACE get_netcdf_variable
+
+    PRIVATE ::  get_netcdf_variable_int, get_netcdf_variable_real
+
  CONTAINS
+
+    SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
+
+       CHARACTER(LEN=PATH), INTENT(IN)         ::  in_file
+       TYPE(nc_var), INTENT(INOUT)             ::  in_var
+       INTEGER(hp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
+
+       INCLUDE 'get_netcdf_variable.inc'
+
+    END SUBROUTINE get_netcdf_variable_int
+
+
+    SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
+
+       CHARACTER(LEN=PATH), INTENT(IN)      ::  in_file
+       TYPE(nc_var), INTENT(INOUT)          ::  in_var
+       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
+
+       INCLUDE 'get_netcdf_variable.inc'
+
+    END SUBROUTINE get_netcdf_variable_real
+
 
     SUBROUTINE netcdf_define_variable(var, ncid)
@@ -59,8 +99,9 @@
 
         CALL check(nf90_def_var(ncid, var % name, NF90_FLOAT,       var % dimids(1:var % ndim), var % varid))
-        CALL check(nf90_put_att(ncid, var % varid, "standard_name", var % standard_name))
         CALL check(nf90_put_att(ncid, var % varid, "long_name",     var % long_name))
         CALL check(nf90_put_att(ncid, var % varid, "units",         var % units))
-        CALL check(nf90_put_att(ncid, var % varid, "lod",           var % lod))
+        IF ( var % lod .GE. 0 )  THEN
+           CALL check(nf90_put_att(ncid, var % varid, "lod",           var % lod))
+        END IF 
         CALL check(nf90_put_att(ncid, var % varid, "source",        var % source))
         CALL check(nf90_put_att(ncid, var % varid, "_FillValue",    NF90_FILL_REAL))
@@ -94,16 +135,10 @@
 !> parameters for the PALM-4U computational grid.
 !------------------------------------------------------------------------------!
-    SUBROUTINE parse_command_line_arguments( start_date, hhl_file,             &
-       soiltyp_file, static_driver_file, input_path, output_file,              &
-       namelist_file, ug, vg, p0, z0, mode )
-
-       CHARACTER(LEN=PATH), INTENT(INOUT)  ::  hhl_file, soiltyp_file,         &
-           static_driver_file, input_path, output_file, namelist_file
-       CHARACTER(LEN=SNAME), INTENT(INOUT) ::  mode
-       REAL(dp), INTENT(INOUT)             ::  ug, vg, p0, z0
-       CHARACTER(LEN=DATE), INTENT(INOUT)  ::  start_date
-
-       CHARACTER(LEN=PATH)     ::  option, arg
-       INTEGER                 ::  arg_count, i
+    SUBROUTINE parse_command_line_arguments( cfg )
+
+       TYPE(inifor_config), INTENT(INOUT) ::  cfg
+
+       CHARACTER(LEN=PATH)                ::  option, arg
+       INTEGER                            ::  arg_count, i
 
        arg_count = COMMAND_ARGUMENT_COUNT()
@@ -111,92 +146,121 @@
 
           ! Every option should have an argument.
-          IF ( MOD(arg_count, 2) .NE. 0 )  THEN
-             message = "Syntax error in command line."
-             CALL abort('parse_command_line_arguments', message)
-          END IF
+          !IF ( MOD(arg_count, 2) .NE. 0 )  THEN
+          !   message = "Syntax error in command line."
+          !   CALL abort('parse_command_line_arguments', message)
+          !END IF
           
           message = "The -clon and -clat command line options are depricated. " // &
              "Please remove them form your inifor command and specify the " // &
              "location of the PALM-4U origin either" // NEW_LINE(' ') // &
-             "   - by setting the namelist parameters 'origin_lon' and 'origin_lat, or'" // NEW_LINE(' ') // &
+             "   - by setting the namelist parameters 'longitude' and 'latitude', or" // NEW_LINE(' ') // &
              "   - by providing a static driver netCDF file via the -static command-line option."
 
-          ! Loop through option/argument pairs.
-          DO i = 1, arg_count, 2
+          i = 1
+          DO WHILE (i .LE. arg_count)
 
              CALL GET_COMMAND_ARGUMENT( i, option )
-             CALL GET_COMMAND_ARGUMENT( i+1, arg )
 
              SELECT CASE( TRIM(option) )
 
-             CASE( '-date' )
-                start_date = TRIM(arg)
+             CASE( '-date', '-d', '--date' )
+                CALL get_option_argument( i, arg )
+                cfg % start_date = TRIM(arg)
 
              ! Elevation of the PALM-4U domain above sea level
-             CASE( '-z0' )
-                READ(arg, *) z0
+             CASE( '-z0', '-z', '--elevation' )
+                CALL get_option_argument( i, arg )
+                READ(arg, *) cfg % z0
 
              ! surface pressure, at z0
-             CASE( '-p0' )
-                READ(arg, *) p0
-
-             ! surface pressure, at z0
-             CASE( '-ug' )
-                READ(arg, *) ug
-
-             ! surface pressure, at z0
-             CASE( '-vg' )
-                READ(arg, *) vg
-
-             ! Domain centre geographical longitude
-             CASE( '-clon' )
+             CASE( '-p0', '-r', '--surface-pressure' )
+                CALL get_option_argument( i, arg )
+                READ(arg, *) cfg % p0
+
+             ! geostrophic wind in x direction
+             CASE( '-ug', '-u', '--geostrophic-u' )
+                CALL get_option_argument( i, arg )
+                READ(arg, *) cfg % ug
+
+             ! geostrophic wind in y direction
+             CASE( '-vg', '-v', '--geostrophic-v' )
+                CALL get_option_argument( i, arg )
+                READ(arg, *) cfg % vg
+
+             ! domain centre geographical longitude and latitude
+             CASE( '-clon', '-clat' )
                 CALL abort('parse_command_line_arguments', message)          
                 !READ(arg, *) lambda_cg
                 !lambda_cg = lambda_cg * TO_RADIANS
-
-             ! Domain centre geographical latitude
-             CASE( '-clat' )
-                CALL abort('parse_command_line_arguments', message)          
                 !READ(arg, *) phi_cg
                 !phi_cg = phi_cg * TO_RADIANS
 
-             CASE( '-path' )
-                 input_path = TRIM(arg)
-
-             CASE( '-hhl' )
-                 hhl_file = TRIM(arg)
-
-             CASE( '-static' )
-                 static_driver_file = TRIM(arg)
-
-             CASE( '-soil' )
-                 soiltyp_file = TRIM(arg)
-
-             CASE( '-o' )
-                output_file = TRIM(arg)
-
-             CASE( '-n' )
-                namelist_file = TRIM(arg)
-
-             ! Initialization mode: 'profile' / 'volume'
-             CASE( '-mode' )
-                mode = TRIM(arg)
-
-                SELECT CASE( TRIM(mode) )
-
-                CASE( 'profile' )
-
-                CASE DEFAULT
-                   message = "Mode '" // TRIM(mode) // "' is not supported. " //&
-                             "Currently, '-mode profile' is the only supported option. " //&
-                             "Select this one or omit the -mode option entirely."
-                   CALL abort( 'parse_command_line_arguments', message ) 
-                END SELECT
+             CASE( '-path', '-p', '--path' )
+                CALL get_option_argument( i, arg )
+                 cfg % input_path = TRIM(arg)
+
+             CASE( '-hhl', '-l', '--hhl-file' )
+                CALL get_option_argument( i, arg )
+                 cfg % hhl_file = TRIM(arg)
+
+             CASE( '-static', '-t', '--static-driver' )
+                CALL get_option_argument( i, arg )
+                 cfg % static_driver_file = TRIM(arg)
+
+             CASE( '-soil', '-s', '--soil-file')
+                CALL get_option_argument( i, arg )
+                 cfg % soiltyp_file = TRIM(arg)
+
+             CASE( '--flow-prefix')
+                CALL get_option_argument( i, arg )
+                 cfg % flow_prefix = TRIM(arg)
+
+             CASE( '--radiation-prefix')
+                CALL get_option_argument( i, arg )
+                 cfg % radiation_prefix = TRIM(arg)
+
+             CASE( '--soil-prefix')
+                CALL get_option_argument( i, arg )
+                 cfg % soil_prefix = TRIM(arg)
+
+             CASE( '--soilmoisture-prefix')
+                CALL get_option_argument( i, arg )
+                 cfg % soilmoisture_prefix = TRIM(arg)
+
+             CASE( '-o', '--output' )
+                CALL get_option_argument( i, arg )
+                cfg % output_file = TRIM(arg)
+
+             CASE( '-n', '--namelist' )
+                CALL get_option_argument( i, arg )
+                cfg % namelist_file = TRIM(arg)
+
+             ! initial condition mode: 'profile' / 'volume'
+             CASE( '-mode', '-i', '--init-mode' )
+                CALL get_option_argument( i, arg )
+                cfg % ic_mode = TRIM(arg)
+
+             ! boundary conditions / forcing mode: 'ideal' / 'real'
+             CASE( '-f', '--forcing-mode' )
+                CALL get_option_argument( i, arg )
+                cfg % bc_mode = TRIM(arg)
+
+             CASE( '--version' )
+                CALL print_version()
+                STOP
+
+             CASE( '--help' )
+                CALL print_version()
+                PRINT *, ""
+                PRINT *, "For a list of command-line options have a look at the README file."
+                STOP
 
              CASE DEFAULT
-                message = "unknown option '" // TRIM(option(2:)) // "'."
+                message = "unknown option '" // TRIM(option) // "'."
                 CALL abort('parse_command_line_arguments', message)
 
              END SELECT
+
+             i = i + 1
 
           END DO
@@ -210,4 +274,75 @@
 
    END SUBROUTINE parse_command_line_arguments
+
+   
+   SUBROUTINE get_option_argument(i, arg)
+      CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
+      INTEGER, INTENT(INOUT)             ::  i
+
+      i = i + 1
+      CALL GET_COMMAND_ARGUMENT(i, arg)
+
+   END SUBROUTINE
+
+
+   SUBROUTINE validate_config(cfg)
+      TYPE(inifor_config), INTENT(IN) ::  cfg
+      LOGICAL                         ::  all_files_present
+
+      all_files_present = .TRUE.
+      all_files_present = all_files_present .AND. file_present(cfg % hhl_file)
+      all_files_present = all_files_present .AND. file_present(cfg % namelist_file)
+      all_files_present = all_files_present .AND. file_present(cfg % output_file)
+      all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file)
+
+      ! Only check optional static driver file name, if it has been given.
+      IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
+         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file)
+      END IF
+
+      IF (.NOT. all_files_present)  THEN
+         message = "INIFOR configuration invalid; some input files are missing."
+         CALL abort( 'validate_config', message ) 
+      END IF
+      
+      
+      SELECT CASE( TRIM(cfg % ic_mode) )
+      CASE( 'profile', 'volume')
+      CASE DEFAULT
+         message = "Initialization mode '" // TRIM(cfg % ic_mode) //&
+                   "' is not supported. " //&
+                   "Please select either 'profile' or 'volume', " //&
+                   "or omit the -i/--init-mode/-mode option entirely, which corresponds "//&
+                   "to the latter."
+         CALL abort( 'validate_config', message ) 
+      END SELECT
+
+
+      SELECT CASE( TRIM(cfg % bc_mode) )
+      CASE( 'real', 'ideal')
+      CASE DEFAULT
+         message = "Forcing mode '" // TRIM(cfg % bc_mode) //& 
+                   "' is not supported. " //&
+                   "Please select either 'real' or 'ideal', " //&
+                   "or omit the -f/--forcing-mode option entirely, which corresponds "//&
+                   "to the latter."
+         CALL abort( 'validate_config', message ) 
+      END SELECT
+
+
+   END SUBROUTINE validate_config
+
+
+   LOGICAL FUNCTION file_present(filename)
+      CHARACTER(LEN=PATH), INTENT(IN) ::  filename
+
+      INQUIRE(FILE=filename, EXIST=file_present)
+
+      IF (.NOT. file_present)  THEN
+         message = "The given file '" // "' does not exist."
+         CALL report('file_present', message)
+      END IF
+
+   END FUNCTION file_present
 
 
@@ -222,15 +357,29 @@
 !> writes the actual data.
 !------------------------------------------------------------------------------!
-    SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid)
+   SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid,                  &
+                                      start_date_string, origin_lon, origin_lat)
 
        TYPE(nc_file), INTENT(INOUT)      ::  output_file
        TYPE(grid_definition), INTENT(IN) ::  palm_grid
-
-       CHARACTER (LEN=SNAME) ::  date
+       CHARACTER (LEN=DATE), INTENT(IN)  ::  start_date_string
+       REAL(dp), INTENT(IN)              ::  origin_lon, origin_lat
+
+       CHARACTER (LEN=8)     ::  date_string
+       CHARACTER (LEN=10)    ::  time_string
+       CHARACTER (LEN=5)     ::  zone_string
+       CHARACTER (LEN=SNAME) ::  history_string
        INTEGER               ::  ncid, nx, ny, nz, nt, dimids(3), dimvarids(3)
        REAL(dp)              ::  z0
 
+       message = "Initializing PALM-4U dynamic driver file '" //               &
+                 TRIM(output_file % name) // "' and setting up dimensions."
+       CALL report('setup_netcdf_dimensions', message)
+
        ! Create the NetCDF file. NF90_CLOBBER selects overwrite mode.
+#if defined( __netcdf4 )
        CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
+#else
+       CALL check(nf90_create(TRIM(output_file % name), NF90_CLOBBER, ncid))
+#endif
 
 !
@@ -238,13 +387,18 @@
 !- Section 1: Write global NetCDF attributes
 !------------------------------------------------------------------------------
-       CALL date_and_time(date)
+       CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
+       history_string =                                                        &
+           'Created on '// date_string      //                                 &
+           ' at '       // time_string(1:2) // ':' // time_string(3:4) //      &
+           ' (UTC'      // zone_string // ')'
+
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title',          'PALM input file for scenario ...'))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution',    'Deutscher Wetterdienst, Offenbach'))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author',         'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
-       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        'Created on '//date))
+       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        TRIM(history_string)))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references',     '--'))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment',        '--'))
-       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     '--'))
-       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     '--'))
+       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
+       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
@@ -267,5 +421,5 @@
           CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
           CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
-          CALL check( nf90_def_dim(ncid, "z", nz+1, dimids(3)) )
+          CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
           output_file % dimids_scl = dimids ! save dimids for later
 
@@ -285,9 +439,9 @@
 
        ! overwrite third dimid with the one of depth
-       CALL check(nf90_def_dim(ncid, "depth", SIZE(palm_grid % depths), dimids(3)) )
+       CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
        output_file % dimids_soil = dimids ! save dimids for later
 
        ! overwrite third dimvarid with the one of depth
-       CALL check(nf90_def_var(ncid, "depth", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
+       CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
        CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
        CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
@@ -301,5 +455,5 @@
           CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
           CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
-          CALL check(nf90_def_dim(ncid, "zw", nz, dimids(3)) )
+          CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
        output_file % dimids_vel = dimids ! save dimids for later
 
@@ -328,5 +482,6 @@
        CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "standard_name", "time"))
        CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "long_name", "time"))
-       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units", "seconds since..."))
+       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units",     &
+                               "seconds since " // start_date_string // " UTC"))
 
        CALL check(nf90_enddef(ncid))
@@ -363,5 +518,5 @@
        INTEGER                              ::  i, ncid
 
-       message = "Initializing PALM-4U forcing file '" // TRIM(filename) // "'."
+       message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
        CALL report('setup_netcdf_variables', message)
 
@@ -374,5 +529,5 @@
 
           IF ( var % to_be_processed )  THEN
-             message = "Defining variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
+             message = "  variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
              CALL report('setup_netcdf_variables', message)
 
@@ -386,5 +541,5 @@
        CALL check(nf90_close(ncid))
 
-       message = "Forcing file '" // TRIM(filename) // "' initialized successfully."
+       message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
        CALL report('setup_netcdf_variables', message)
 
@@ -447,5 +602,5 @@
 
           input_var => group % in_var_list(1)
-          buffer(buf_id) % array = get_netcdf_variable( input_file, input_var )  
+          CALL get_netcdf_variable(input_file, input_var, buffer(buf_id) % array)
           CALL report('read_input_variables', "Read accumulated " // TRIM(group % in_var_list(1) % name)) 
 
@@ -472,5 +627,5 @@
              END IF
 
-             buffer(ivar) % array = get_netcdf_variable( input_file, input_var )  
+             CALL get_netcdf_variable(input_file, input_var, buffer(ivar) % array)
 
              IF ( input_var % is_upside_down )  CALL reverse(buffer(ivar) % array)
@@ -545,4 +700,5 @@
 
           CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
+          CALL check(nf90_close(ncid))
 
        ELSE
@@ -555,90 +711,4 @@
 
     END FUNCTION get_netcdf_attribute
-
-
-
-    FUNCTION get_netcdf_variable(in_file, in_var) RESULT(buffer)
-
-       CHARACTER(LEN=PATH), INTENT(IN)      ::  in_file
-       TYPE(nc_var), INTENT(INOUT)          ::  in_var
-       REAL(dp), ALLOCATABLE                ::  buffer(:,:,:)
-       INTEGER                              ::  i, ncid, start(3)
-
-
-       ! Read in_var NetCDF attributes
-       IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
-            nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR )  THEN
-
-          CALL check(nf90_get_att(ncid, in_var % varid, "long_name", in_var % long_name))
-          CALL check(nf90_get_att(ncid, in_var % varid, "units", in_var % units))
-
-          ! Read in_var NetCDF dimensions
-          CALL check(nf90_inquire_variable( ncid, in_var % varid,              &
-                                            ndims  = in_var % ndim,            &
-                                            dimids = in_var % dimids ))
-                                        
-          DO i = 1, in_var % ndim
-             CALL check(nf90_inquire_dimension( ncid, in_var % dimids(i),      &
-                                                name = in_var % dimname(i),    &
-                                                len  = in_var % dimlen(i) ))
-          END DO
-
-          start = (/ 1, 1, 1 /)
-          IF ( TRIM(in_var % name) .EQ. 'T_SO' )  THEN
-             ! Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
-             in_var % dimlen(3) = in_var % dimlen(3) - 1
-
-             ! Start reading from second level, e.g. depth = 0.005 instead of 0.0
-             start(3) = 2
-          END IF
-
-          SELECT CASE(in_var % ndim)
-
-          CASE (2)
-
-             ALLOCATE( buffer( in_var % dimlen(1),                             &
-                               in_var % dimlen(2),                             &
-                               1 ) )
-
-          CASE (3)
-
-             ALLOCATE( buffer( in_var % dimlen(1),                             &
-                               in_var % dimlen(2),                             &
-                               in_var % dimlen(3) ) )
-          CASE (4)
-
-             ALLOCATE( buffer( in_var % dimlen(1),                             &
-                               in_var % dimlen(2),                             &
-                               in_var % dimlen(3) ) )
-          CASE DEFAULT
-
-             message = "Failed reading NetCDF variable " //                    &
-                TRIM(in_var % name) // " with " // TRIM(str(in_var%ndim)) //   &
-                " dimensions because only two- and and three-dimensional" //   &
-                " variables are supported."
-             CALL abort('get_netcdf_variable', message) 
-
-          END SELECT
- CALL run_control('time', 'alloc')
-          
-          ! TODO: Check for matching dimensions of buffer and var
-          CALL check(nf90_get_var( ncid, in_var % varid, buffer,               &
-                                   start = start,                              &
-                                   count = in_var % dimlen(1:3) ) )
-
- CALL run_control('time', 'read')
-       ELSE
-
-          message = "Failed to read '" // TRIM(in_var % name) // &
-             "' from file '" // TRIM(in_file) // "'."
-          CALL report('get_netcdf_variable', message)
-
-       END IF
-
-       CALL check(nf90_close(ncid))
-
- CALL run_control('time', 'read')
-
-    END FUNCTION get_netcdf_variable
 
 
@@ -657,9 +727,6 @@
 
        ! Skip time dimension for output
-       IF ( var_is_time_dependent )  THEN
-          ndim = var % ndim - 1
-       ELSE 
-          ndim = var % ndim
-       END IF
+       ndim = var % ndim
+       IF ( var_is_time_dependent )  ndim = var % ndim - 1
 
        start(:)      = (/1,1,1,1/)
@@ -733,9 +800,14 @@
                                    start=start(1:ndim+1) ) )
 
-       CASE ( 'profile' )
+       CASE ( 'constant scalar profile' )
 
           CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
                                    start=start(1:ndim+1),                      &
                                    count=count(1:ndim) ) )
+
+       CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
+
+           message = "Doing nothing in terms of writing large-scale forings."
+           CALL report('update_output', message)
 
        CASE DEFAULT
Index: palm/trunk/UTIL/inifor/src/transform.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/transform.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/transform.f90	(revision 3182)
@@ -21,4 +21,7 @@
 ! Current revisions:
 ! -----------------
+! Introduced new PALM grid stretching
+! Removed unnecessary subroutine parameters
+! Renamed kcur to k_intermediate
 ! 
 ! 
@@ -80,15 +83,13 @@
        TYPE(grid_definition), INTENT(IN) ::  outgrid
        REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
-       REAL(dp), INTENT(OUT)             ::  out_arr(0:,0:,0:)
-
-       INTEGER :: i, j, k, l, nx, ny, nz
-
-       nx = UBOUND(out_arr, 1)
-       ny = UBOUND(out_arr, 2)
+       REAL(dp), INTENT(OUT)             ::  out_arr(0:,0:,:)
+
+       INTEGER :: i, j, k, l, nz
+
        nz = UBOUND(out_arr, 3)
 
-       DO j = 0, ny
-       DO i = 0, nx
-       DO k = nz, 0, -1
+       DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2)
+       DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1)
+       DO k = nz, LBOUND(out_arr, 3), -1
 
           ! TODO: Remove IF clause and extrapolate based on a critical vertical
@@ -101,6 +102,6 @@
              out_arr(i,j,k) = 0.0_dp
              DO l = 1, 2
-                out_arr(i,j,k) = out_arr(i,j,k) +                                 &
-                    outgrid % w_verti(i,j,k,l) *                                  &
+                out_arr(i,j,k) = out_arr(i,j,k) +                              &
+                    outgrid % w_verti(i,j,k,l) *                               &
                     in_arr(i,j,outgrid % kk(i,j,k, l) )
              END DO
@@ -139,7 +140,7 @@
     ! I index 0-based for the indices of the outvar to be consistent with the
     ! outgrid indices and interpolation weights.
-       TYPE(grid_definition), INTENT(IN) ::  outgrid
-       REAL(dp), INTENT(IN)              ::  invar(0:,0:,0:)
-       REAL(dp), INTENT(OUT)             ::  outvar(0:,0:,0:)
+       TYPE(grid_definition), INTENT(IN)  ::  outgrid
+       REAL(dp), INTENT(IN)               ::  invar(0:,0:,0:)
+       REAL(dp), INTENT(OUT)              ::  outvar(0:,0:,0:)
        TYPE(nc_var), INTENT(IN), OPTIONAL ::  ncvar
 
@@ -413,4 +414,5 @@
 
     END SUBROUTINE rotate_to_cosmo
+       
 
 
@@ -427,6 +429,6 @@
 !>                     -------------
 !>           jj, lat
-!>              ^      j
-!>              |       \          i
+!>              ^        j
+!>              |         \          i
 !>  jj(i,j,2/3) + ... 2 ---\--------/------ 3
 !>              |     | ^   \      /        |
@@ -459,10 +461,11 @@
 !> 
 !------------------------------------------------------------------------------!
-    SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, cosmo_dxi,     &
-       cosmo_dyi, palm_clat, palm_clon, palm_ii, palm_jj)
+    SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon,                &
+                                          palm_clat, palm_clon,                &
+                                          palm_ii, palm_jj)
 
        REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
        REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
-       REAL(dp), INTENT(IN)                       ::  cosmo_dxi, cosmo_dyi
+       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
        INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT)  ::  palm_ii, palm_jj
 
@@ -472,4 +475,6 @@
        lon0 = cosmo_lon(0)
        lat0 = cosmo_lat(0)
+       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
+       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
 
        DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
@@ -508,5 +513,5 @@
        TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
 
-       INTEGER  ::  i, j, k, nx, ny, nz, nlev, kcur
+       INTEGER  ::  i, j, k, nx, ny, nz, nlev, k_intermediate
        LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
                     point_is_in_current_cell
@@ -523,5 +528,5 @@
        DO j = 0, ny
 
-          kcur = 0
+          k_intermediate = 0
 
           column_base = palm_intermediate % h(i,j,0)
@@ -532,11 +537,11 @@
           ! cell, or above column_top. Keep increasing current cell index until
           ! the current cell overlaps with the current_height.
-          DO k = 0, nz
+          DO k = 1, nz
 
              ! Memorize the top and bottom boundaries of the coarse cell and the
              ! current height within it
              current_height = palm_grid % z(k) + palm_grid % z0
-             h_top    = palm_intermediate % h(i,j,kcur+1)
-             h_bottom = palm_intermediate % h(i,j,kcur)
+             h_top    = palm_intermediate % h(i,j,k_intermediate+1)
+             h_bottom = palm_intermediate % h(i,j,k_intermediate)
 
              point_is_above_grid = (current_height > column_top) !22000m, very unlikely
@@ -556,4 +561,7 @@
                 palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
 
+                message = "PALM-4U grid extends above COSMO-DE model top."
+                CALL abort('find_vertical_neighbours_and_weights', message)
+
              ELSE IF (point_is_below_grid)  THEN
 
@@ -564,9 +572,9 @@
                 ! cycle through intermediate levels until current
                 ! intermediate-grid cell overlaps with current_height
-                DO WHILE (.NOT. point_is_in_current_cell .AND. kcur <= nlev-1)
-                   kcur = kcur + 1
-
-                   h_top    = palm_intermediate % h(i,j,kcur+1)
-                   h_bottom = palm_intermediate % h(i,j,kcur)
+                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
+                   k_intermediate = k_intermediate + 1
+
+                   h_top    = palm_intermediate % h(i,j,k_intermediate+1)
+                   h_bottom = palm_intermediate % h(i,j,k_intermediate)
                    point_is_in_current_cell = (                                &
                       current_height >= h_bottom .AND.                         &
@@ -575,13 +583,14 @@
                 END DO
 
-                ! kcur = 48 indicates the last section (indices 48 and 49), i.e.
-                ! kcur = 49 is not the beginning of a valid cell.
-                IF (kcur > nlev-1)  THEN
-                   message = "Index " // TRIM(str(kcur)) // " is above intermediate grid range."
+                ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
+                ! k_intermediate = 49 is not the beginning of a valid cell.
+                IF (k_intermediate > nlev-1)  THEN
+                   message = "Index " // TRIM(str(k_intermediate)) //          &
+                             " is above intermediate grid range."
                    CALL abort('find_vertical_neighbours', message)
                 END IF
    
-                palm_grid % kk(i,j,k,1) = kcur
-                palm_grid % kk(i,j,k,2) = kcur + 1
+                palm_grid % kk(i,j,k,1) = k_intermediate
+                palm_grid % kk(i,j,k,2) = k_intermediate + 1
 
                 ! copmute vertical weights
@@ -643,8 +652,8 @@
 !          
     SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
-       cosmo_dxi, cosmo_dyi, palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
+       palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
        
        REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
-       REAL(dp), INTENT(IN)                       ::  cosmo_dxi, cosmo_dyi
+       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
        REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
        INTEGER, DIMENSION(0:,0:,1:), INTENT(IN)   ::  palm_ii, palm_jj
@@ -654,4 +663,7 @@
        REAL(dp) ::  wl, wp
        INTEGER  ::  i, j
+
+       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
+       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
 
        DO j = 0, UBOUND(palm_clon, 2)
Index: palm/trunk/UTIL/inifor/src/types.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/types.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/types.f90	(revision 3182)
@@ -21,4 +21,9 @@
 ! Current revisions:
 ! -----------------
+! Introduced new PALM grid stretching:
+! - Converted vertical grid_definition coordinte variables to pointers
+! Improved command line interface:
+! - Moved INIFOR configuration into a new derived data type
+! Removed unnecessary variables 
 ! 
 ! 
@@ -41,9 +46,34 @@
  
  USE defs,                                                                     &
-    ONLY:  dp, PATH, SNAME, LNAME
+    ONLY:  dp, DATE, PATH, SNAME, LNAME
  USE netcdf,                                                                   &
     ONLY:  NF90_MAX_VAR_DIMS, NF90_MAX_NAME
 
  IMPLICIT NONE
+
+ TYPE inifor_config
+    CHARACTER(LEN=DATE)  ::  start_date           !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation
+
+    CHARACTER(LEN=PATH)  ::  input_path           !< Path to the input data file directory
+    CHARACTER(LEN=PATH)  ::  hhl_file             !< Path to the file containing the COSMO-DE HHL variable (height of half layers, i.e. vertical cell faces)
+    CHARACTER(LEN=PATH)  ::  namelist_file        !< Path to the PALM-4U namelist file
+    CHARACTER(LEN=PATH)  ::  output_file          !< Path to the INIFOR output file (i.e. PALM-4U dynamic driver')
+    CHARACTER(LEN=PATH)  ::  soiltyp_file         !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
+    CHARACTER(LEN=PATH)  ::  static_driver_file   !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)
+
+    CHARACTER(LEN=SNAME) ::  flow_prefix          !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  soil_prefix          !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses
+
+    CHARACTER(LEN=SNAME) ::  bc_mode
+    CHARACTER(LEN=SNAME) ::  ic_mode
+    CHARACTER(LEN=SNAME) ::  rotation_method
+
+    REAL(dp)             ::  p0
+    REAL(dp)             ::  ug
+    REAL(dp)             ::  vg
+    REAL(dp)             ::  z0 !< Elevation of the PALM-4U domain above sea level [m]
+ END TYPE inifor_config
 
  TYPE grid_definition
@@ -55,13 +85,6 @@
     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.
     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.
-    REAL(dp)              ::  dx            !< grid spacing in the first dimension [m]
-    REAL(dp)              ::  dy            !< grid spacing in the second dimension [m]
-    REAL(dp)              ::  dz            !< grid spacing in the third dimension [m]
-    REAL(dp)              ::  dxi           !< inverse grid spacing in the first dimension [m^-1]
-    REAL(dp)              ::  dyi           !< inverse grid spacing in the second dimension [m^-1]
-    REAL(dp)              ::  dzi           !< inverse grid spacing in the third dimension [m^-1]
     REAL(dp)              ::  lx            !< domain length in the first dimension [m]
     REAL(dp)              ::  ly            !< domain length in the second dimension [m]
-    REAL(dp)              ::  lz            !< domain length in the third dimension [m]
     REAL(dp)              ::  x0            !< x coordinate of PALM-4U domain projection centre, i.e. location of zero distortion
     REAL(dp)              ::  y0            !< y coordinate of PALM-4U domain projection centre, i.e. location of zwro distortion
@@ -69,5 +92,5 @@
     REAL(dp), ALLOCATABLE ::  x(:)          !< coordinates of cell centers in x direction [m]
     REAL(dp), ALLOCATABLE ::  y(:)          !< coordinates of cell centers in y direction [m]
-    REAL(dp), ALLOCATABLE ::  z(:)          !< coordinates of cell centers in z direction [m]
+    REAL(dp), POINTER     ::  z(:)          !< coordinates of cell centers in z direction [m]
     REAL(dp), ALLOCATABLE ::  h(:,:,:)      !< heights grid point for intermediate grids [m]
     REAL(dp), POINTER     ::  hhl(:,:,:)    !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from 
@@ -76,5 +99,5 @@
     REAL(dp), ALLOCATABLE ::  xu(:)         !< coordinates of cell faces in x direction [m]
     REAL(dp), ALLOCATABLE ::  yv(:)         !< coordinates of cell faces in y direction [m]
-    REAL(dp), ALLOCATABLE ::  zw(:)         !< coordinates of cell faces in z direction [m]
+    REAL(dp), POINTER     ::  zw(:)         !< coordinates of cell faces in z direction [m]
     REAL(dp), ALLOCATABLE ::  lat(:)        !< rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]
     REAL(dp), ALLOCATABLE ::  lon(:)        !< rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]
@@ -89,5 +112,5 @@
     REAL(dp), ALLOCATABLE ::  w_horiz(:,:,:)   !< weights for bilinear horizontal interpolation
     REAL(dp), ALLOCATABLE ::  w_verti(:,:,:,:) !< weights for linear vertical interpolation
- END TYPE
+ END TYPE grid_definition
 
 
@@ -103,5 +126,5 @@
     INTEGER               ::  dimvarids_soil(3)!< NetCDF IDs of the grid coordinates for soil points x, y, depth 
     REAL(dp), POINTER     ::  time(:)       ! vector of output time steps
- END TYPE
+ END TYPE nc_file
 
 
@@ -123,7 +146,7 @@
     CHARACTER(LEN=SNAME)                  ::  kind                      !< Kind of grid
     CHARACTER(LEN=SNAME)                  ::  task                      !< Processing task that generates this variable, e.g. 'interpolate_2d' or 'average profile'
-    LOGICAL                               ::  to_be_processed = .FALSE. !< Inifor flag indicating whether variable shall be processed
-    LOGICAL                               ::  is_read = .FALSE.         !< Inifor flag indicating whether variable has been read
-    LOGICAL                               ::  is_upside_down  = .FALSE. !< Inifor flag indicating whether variable shall be processed
+    LOGICAL                               ::  to_be_processed = .FALSE. !< INIFOR flag indicating whether variable shall be processed
+    LOGICAL                               ::  is_read = .FALSE.         !< INIFOR flag indicating whether variable has been read
+    LOGICAL                               ::  is_upside_down  = .FALSE. !< INIFOR flag indicating whether variable shall be processed
     TYPE(grid_definition), POINTER        ::  grid                      !< Pointer to the corresponding output grid
     TYPE(grid_definition), POINTER        ::  intermediate_grid         !< Pointer to the corresponding intermediate grid
Index: palm/trunk/UTIL/inifor/src/util.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/util.f90	(revision 2718)
+++ palm/trunk/UTIL/inifor/src/util.f90	(revision 3182)
@@ -21,4 +21,5 @@
 ! Current revisions:
 ! -----------------
+! Improved real-to-string conversion
 ! 
 ! 
@@ -43,5 +44,5 @@
         ONLY :  C_CHAR, C_INT, C_PTR, C_SIZE_T
     USE defs,                                                                  &
-        ONLY :  dp, PI, DATE
+        ONLY :  dp, PI, DATE, SNAME
 
     IMPLICIT NONE
@@ -279,5 +280,5 @@
     ! Convert a real number to a string in scientific notation
     ! showing four significant digits.
-    CHARACTER(LEN=11) FUNCTION real_to_str(val, format)
+    CHARACTER(LEN=SNAME) FUNCTION real_to_str(val, format)
 
         REAL(dp), INTENT(IN)                   ::  val
@@ -285,9 +286,9 @@
 
         IF (PRESENT(format))  THEN
-           WRITE(real_to_str, TRIM(format)) val
+           WRITE(real_to_str, format) val
         ELSE
            WRITE(real_to_str, '(E11.4)') val
-           real_to_str = ADJUSTL(real_to_str)
         END IF
+        real_to_str = ADJUSTL(real_to_str)
 
     END FUNCTION real_to_str
