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')
 
