Index: /palm/trunk/UTIL/inifor/Makefile
===================================================================
--- /palm/trunk/UTIL/inifor/Makefile	(revision 3394)
+++ /palm/trunk/UTIL/inifor/Makefile	(revision 3395)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Updated build order
+# 
+# 3183 2018-07-27 14:25:55Z suehring
 # Added __netcdf4 preprocessor flag
 # 
@@ -47,5 +50,5 @@
 TEST_PATH = $(PROJECT_PATH)/tests
 
-MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
+MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
 SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90
 OBJECTS = $(SOURCES:%.f90=%.o)
Index: /palm/trunk/UTIL/inifor/Makefile.gnu
===================================================================
--- /palm/trunk/UTIL/inifor/Makefile.gnu	(revision 3394)
+++ /palm/trunk/UTIL/inifor/Makefile.gnu	(revision 3395)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Updated build order
+# 
+# 3183 2018-07-27 14:25:55Z suehring
 # Added __netcdf4 preprocessor flag
 # Corrected compilation order
@@ -44,9 +47,9 @@
 PROJECT = inifor
 PROJECT_PATH = .
-BIN_PATH  = $(PROJECT_PATH)/bin
+BIN_PATH  = $(PROJECT_PATH)/../../SCRIPTS
 SRC_PATH  = $(PROJECT_PATH)/src
 TEST_PATH = $(PROJECT_PATH)/tests
 
-MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
+MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
 SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90
 OBJECTS = $(SOURCES:%.f90=%.o)
Index: /palm/trunk/UTIL/inifor/Makefile.ifort
===================================================================
--- /palm/trunk/UTIL/inifor/Makefile.ifort	(revision 3394)
+++ /palm/trunk/UTIL/inifor/Makefile.ifort	(revision 3395)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Updated build order
+# 
+# 3183 2018-07-27 14:25:55Z suehring
 # Added __netcdf4 preprocessor flag
 # Corrected compilation order
@@ -44,9 +47,9 @@
 PROJECT = inifor
 PROJECT_PATH = .
-BIN_PATH  = $(PROJECT_PATH)/bin
+BIN_PATH  = $(PROJECT_PATH)/../../SCRIPTS
 SRC_PATH  = $(PROJECT_PATH)/src
 TEST_PATH = $(PROJECT_PATH)/tests
 
-MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
+MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
 SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90
 OBJECTS = $(SOURCES:%.f90=%.o)
Index: /palm/trunk/UTIL/inifor/README
===================================================================
--- /palm/trunk/UTIL/inifor/README	(revision 3394)
+++ /palm/trunk/UTIL/inifor/README	(revision 3395)
@@ -1,3 +1,3 @@
-# INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.3.0)
+# INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.4.0)
 
 INIFOR provides the meteorological fields required to initialize and drive the
@@ -14,11 +14,55 @@
 
 
+## MINIMAL USAGE EXAMPLE
+
+1. Customize `./namelist` (number of grid points and spacings, end_time)
+2. Run `current_version/trunk/SCRIPTS/inifor --path <scenario path> --date <YYYYMMDD>`
+
+
 ## USAGE
 
-1. Customize `./namelist` (number or grid points and spacings, end_time)
-2. Run `current_version/trunk/SCRIPTS/inifor -path <scenario path> -date <YYYYMMDD>`
+After compilation, the `inifor` binary resides in the `$PALM_BIN` path, i.e. in
+`<PALM path>/current_version/trunk/SCRIPTS/`.
+
+In order to run, INIFOR requires three kinds of inputs:
+
+1. hourly COSMO model output,
+2. a steering namelist file, and
+3. command-line options.
+
+In addition, a static driver file may be supplied (see `--static-driver` option
+below) in order to pass the coordinates of the PALM origin to INIFOR. If no
+static driver is passed to INIFOR, origin coordinates are read from the namelist
+file.
+
+A typical `inifor` call looks like this:
+
+    inifor --path /data/evaluation/20170729 --date 2017073006 \
+           --init-mode profile -n namelist -o dynamic_driver.nc \
+           --elevation 110 --input-prefix lff0
+
+### INPUT DATA: COSMO MODEL OUTPUT
+
+INIFOR processes COSMO model output which it requires to be stored in a set of
+netCDF files located in a user-specified path (see `--path` option). These are:
+
+- hhl.nc: This file provides the COSMO numerical grid. (*hhl* abbreviates
+  *height of half layers*, i.e. the heights of the vertical cell boundaries.)
+- soil.nc: This file provides the COSMO soil map and is used to destinguish land
+  from water cells.
+- `<prefix>YYYYMMDDHH-<suffix>.nc`: Each of these files contains COSMO model
+  of one time step at the given time in UTC (Y..year, M..month, D..day, H..hour).
+    - The `<prefix>` distinguishes different DWD products, for instance COSMO
+      analyses (`laf`) or forecasts (`lff`).
+    - The `<suffix>` distinguishes four kinds of COSMO model output data, namely
+        - `flow` (atmospheric fields)
+        - `rad` (radiation)
+        - `soil` (soil moisture and temperature)
+        - `soilmoisture` (precipitation and evaporation)
+    - For example, laf2016010100-flow.nc contains the atmospheric fields of the
+      COSMO analysis of Januray 1st, 2016 for 0:00 UTC.
 
 
-## AVAILABLE NAMELIST PARAMETERS
+### AVAILABLE NAMELIST PARAMETERS
 
 INIFOR mirrors a subset of the PALM-4U's Fortran namelists `inipar` and `d3par`
@@ -26,5 +70,5 @@
 
 
-### inipar
+#### inipar
 
     nx, ny, nz - number of PALM-4U grid points in x, y, and z direction
@@ -34,11 +78,11 @@
     dz_max - maximum vertical grid spacing [m]
     dz_stretch_level_start(9) - array of height levels above which the grid is
-    to be stretched vertically [m]
+        to be stretched vertically [m]
     dz_stretch_level_end(9) - array of height levels until which the grid is to
-    be stretched vertically [m]
+        be stretched vertically [m]
     longitude, latitude - geographical coordinates of the PALM-4U origin [deg]
 
 
-### d3par
+#### d3par
 
     end_time - PALM-4U simulation time. INIFOR will produce hourly forcing data
@@ -47,5 +91,5 @@
 
 
-### EXAMPLE NAMELIST FILE
+#### EXAMPLE NAMELIST FILE
 
     &inipar nx = 4679, ny = 3939, nz = 360
@@ -59,5 +103,13 @@
 
 
-## AVAILABLE COMMAND-LINE PARAMETERS
+### AVAILABLE COMMAND-LINE PARAMETERS
+    --averaging-mode <mode>:
+        Selects how averaged quantities (large-scale forcing terms) are computed.
+        INIFOR supports averaging along input model levels ('level') and along
+        constant heights ('height'). Default: level
+
+    -a, --averaging-angle <angle>:
+        Width of the averaging box in longitudal and latitudal direction in the
+        source coordinate system (COSMO rotated-pole) [deg]. Default: 2.0
 
     -d, --date, -date <date>:
@@ -65,5 +117,5 @@
         hours (HH) are given, INIFOR assumes that the simulation starts at O UTC
         on that day. Default: 20130721
-        
+
     -i, --init-mode, -mode <mode>:
         Set the PALM-4U initialization mode. INIFOR can provide initial conditions
@@ -72,5 +124,5 @@
 
     -l, --hhl-file, -hhl <netCDF file>:
-        Location of the netCDF file containing the vertical COSMO-DE grid levels,
+        Location of the netCDF file containing the vertical COSMO grid levels,
         specifically the heights of half levels (hhl, i.e. vertical cell faces).
         Default: <scenario path>/hhl.nc
@@ -83,5 +135,5 @@
     -o, --output <output file>:
         Name of the INIFOR output file, i.e. the PALM-4U dynamic driver.
-        Default: ./palm-4u-input.nc
+        Default: ./dynamic_driver.nc
 
     -p, --path, -path <scenario path>:
@@ -89,23 +141,22 @@
 
     -r, --surface-pressure, -p0 <pressure>:
-        Surface pressure at z=0 in the PALM-4U domain [Pa].
-        Default: 1e5
+        Manually set the pressure at z=0 in the PALM-4U domain [Pa]. If not
+        given, the surface pressure is computed from the COSMO pressure field.
 
     -s, --soil-file, -soil <netCDF file>:
-        Location of the netCDF file containing the COSMO-DE soil type map.
+        Location of the netCDF file containing the COSMO soil type map.
         Default: <scenario path>/soil.nc
 
     -t, --static-driver, -static <netCDF file>:
-        Location of the netCDF file containing the static driver file for the case
+        Location of the netCDF file containing the static driver for the case
         to be simulated with PALM-4U. Optional parameter. Default: None
 
     -u, --geostrophic-u, -ug <velocity>:
-        Specifies the geostrophic wind in x direction [m/s]. Default: 0
+        Manually specify the geostrophic wind in x direction [m/s]. If not
+        given, the geostrophic wind is computed from the COSMO pressure field.
 
     -v, --geostrophic-v, -vg <velocity>:
-        Specifies the geostrophic wind in y direction [m/s]. Default: 0
-
-    --version:
-        Output version number and exit.
+        Manually specify the geostrophic wind in y direction [m/s]. If not
+        given, the geostrophic wind is computed from the COSMO pressure field.
 
     -z, --elevation, -z0 <height>: Specifies the elevation of the PALM-4U domain
@@ -113,5 +164,9 @@
 
 
-## ADDITIONAL COMMAND-LINE OPTIONS
+#### ADDITIONAL COMMAND-LINE OPTIONS
+
+    --input-prefix <prefix>:
+        Set the file prefixes for all input files. Individual prefixes can be
+        overwritten with the options below. Default: laf
 
     --flow-prefix <prefix>:
@@ -126,2 +181,9 @@
     --soilmoisture-prefix <prefix>:
         Set the file prefix of soil moisture input files. Default: laf
+
+    --debug:
+        Enable debugging messages.
+
+    --version:
+        Output version number and exit.
+
Index: /palm/trunk/UTIL/inifor/src/control.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/control.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/control.f90	(revision 3395)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Suppress debugging messages unless --debug option is given
+! 
+!
+! 3183 2018-07-27 14:25:55Z suehring
 ! Added version and copyright output
 ! 
@@ -57,12 +61,13 @@
  CONTAINS
 
-    SUBROUTINE report(routine, message)
-
-       CHARACTER(LEN=*), INTENT(IN) ::  routine
-       CHARACTER(LEN=*), INTENT(IN) ::  message
-       INTEGER                      ::  u
-       LOGICAL, SAVE                ::  is_first_run = .TRUE.
-
-       PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
+    SUBROUTINE report(routine, message, debug)
+
+       CHARACTER(LEN=*), INTENT(IN)  ::  routine
+       CHARACTER(LEN=*), INTENT(IN)  ::  message
+       LOGICAL, OPTIONAL, INTENT(IN) ::  debug
+       INTEGER                       ::  u
+       LOGICAL, SAVE                 ::  is_first_run = .TRUE.
+       LOGICAL                       ::  suppress_message
+
 
        IF (is_first_run)  THEN
@@ -73,5 +78,14 @@
        END IF
           
-       WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
+
+       suppress_message = .FALSE.
+       IF (PRESENT(debug))  THEN
+          IF (.NOT. debug)  suppress_message = .TRUE.
+       END IF
+
+       IF (.NOT. suppress_message)  THEN
+          PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
+          WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
+       END IF
 
        CLOSE(u)
Index: /palm/trunk/UTIL/inifor/src/defs.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/defs.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/defs.f90	(revision 3395)
@@ -26,4 +26,9 @@
 ! -----------------
 ! $Id$
+! New parameters for computation of geostrophic winds
+! Bumped INIFOR version number
+! 
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Updated defaults for soil extrapolation steps and nudging time-scale
 ! Improved handling of the start date string
@@ -73,6 +78,8 @@
  REAL(dp), PARAMETER ::  G    = 9.80665_dp                    !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s]
  REAL(dp), PARAMETER ::  RHO_L = 1e3_dp                       !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3], in [kg/m^3]
+ REAL(dp), PARAMETER ::  HECTO = 100_dp                       !< unit conversion factor from hPa to Pa
 
  ! PALM-4U parameters
+ REAL(dp), PARAMETER ::  OMEGA   = 7.29e-5_dp                 !< angular velocity of Earth's rotation [s^-1]
  REAL(dp), PARAMETER ::  P_REF   = 1e5_dp                     !< Reference pressure for potential temperature [Pa]
  REAL(dp), PARAMETER ::  RD_PALM = 287.0_dp                   !< specific gas constant of dry air, used in computation of PALM-4U's potential temperature [J/kg/K]
@@ -83,5 +90,5 @@
  INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
  REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
- CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.3.0'            !< INIFOR version number
+ CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.0'            !< INIFOR version number
  CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2018 Leibniz Universitaet Hannover' // &
      NEW_LINE(' ') // ' Copyright 2017-2018 Deutscher Wetterdienst Offenbach' !< Copyright notice
Index: lm/trunk/UTIL/inifor/src/get_netcdf_variable.inc
===================================================================
--- /palm/trunk/UTIL/inifor/src/get_netcdf_variable.inc	(revision 3394)
+++ 	(revision )
@@ -1,118 +1,0 @@
-!> @file src/get_netcdf_variable.inc
-!------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017-2018 Leibniz Universitaet Hannover
-! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
-!------------------------------------------------------------------------------!
-!
-! Current revisions:
-! -----------------
-! 
-! 
-! Former revisions:
-! -----------------
-! $Id$
-! Initial revision
-! 
-! 
-! 
-!
-! 
-!
-! Authors:
-! --------
-! @author Eckhard Kadasch
-!
-! Description:
-! ------------
-!> This file provides the implementation of the io.f90 subroutines
-!> get_netcdf_variable_int() and get_netcdf_variable_real(). Both routines
-!> differ only in the type of the netCDF variable to be read and otherwise share
-!> their functionality. In order to avoid duplicating code, both routines source
-!> the implementation below using Fortran's INCLUDE statement.
-!> 
-!------------------------------------------------------------------------------!
- INTEGER ::  i, ncid, start(3)
-
- ! Read in_var NetCDF attributes
- IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
-      nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR )  THEN
-
-    CALL check(nf90_get_att(ncid, in_var % varid, "long_name", in_var % long_name))
-    CALL check(nf90_get_att(ncid, in_var % varid, "units", in_var % units))
-
-    ! Read in_var NetCDF dimensions
-    CALL check(nf90_inquire_variable( ncid, in_var % varid,              &
-                                      ndims  = in_var % ndim,            &
-                                      dimids = in_var % dimids ))
-                                  
-    DO i = 1, in_var % ndim
-       CALL check(nf90_inquire_dimension( ncid, in_var % dimids(i),      &
-                                          name = in_var % dimname(i),    &
-                                          len  = in_var % dimlen(i) ))
-    END DO
-
-    start = (/ 1, 1, 1 /)
-    IF ( TRIM(in_var % name) .EQ. 'T_SO' )  THEN
-       ! Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
-       in_var % dimlen(3) = in_var % dimlen(3) - 1
-
-       ! Start reading from second level, e.g. depth = 0.005 instead of 0.0
-       start(3) = 2
-    END IF
-
-    SELECT CASE(in_var % ndim)
-
-    CASE (2)
-
-       ALLOCATE( buffer( in_var % dimlen(1),                             &
-                         in_var % dimlen(2),                             &
-                         1 ) )
-
-    CASE (3, 4)
-
-       ALLOCATE( buffer( in_var % dimlen(1),                             &
-                         in_var % dimlen(2),                             &
-                         in_var % dimlen(3) ) )
-
-    CASE DEFAULT
-
-       message = "Failed reading NetCDF variable " //                    &
-          TRIM(in_var % name) // " with " // TRIM(str(in_var%ndim)) //   &
-          " dimensions because only two- and and three-dimensional" //   &
-          " variables are supported."
-       CALL abort('get_netcdf_variable', message) 
-
-    END SELECT
- CALL run_control('time', 'alloc')
-          
-    ! TODO: Check for matching dimensions of buffer and var
-    CALL check(nf90_get_var( ncid, in_var % varid, buffer,               &
-                             start = start,                              &
-                             count = in_var % dimlen(1:3) ) )
-
- CALL run_control('time', 'read')
- ELSE
-
-    message = "Failed to read '" // TRIM(in_var % name) // &
-       "' from file '" // TRIM(in_file) // "'."
-    CALL report('get_netcdf_variable', message)
-
- END IF
-
- CALL check(nf90_close(ncid))
-
- CALL run_control('time', 'read')
Index: /palm/trunk/UTIL/inifor/src/grid.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/grid.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/grid.f90	(revision 3395)
@@ -26,4 +26,14 @@
 ! -----------------
 ! $Id$
+! Added computation of geostrophic winds form COSMO pressure fields
+! Introduced averaging grids and internal 'output' variables for computation of
+!     geostrophic and large-scale forcings
+! Merged qv, uv and vg output variables and 'temperature' IO group into new
+!     'thermodynamics' IP group
+! Cleaned up terminal output, show some messages only with --debug option
+! Improved variable naming and readability
+! 
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Introduced new PALM grid stretching
 ! Updated variable names and metadata for PIDS v1.9 compatibility
@@ -61,5 +71,6 @@
         ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp,     &
                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
+               BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM,     &
+               RHO_L, OMEGA, HECTO
     USE io,                                                                    &
         ONLY:  get_netcdf_variable, get_netcdf_attribute,                      &
@@ -68,9 +79,11 @@
         ONLY:  NF90_MAX_NAME, NF90_MAX_VAR_DIMS
     USE transform,                                                             &
-        ONLY:  rotate_to_cosmo, find_horizontal_neighbours,                    &
+        ONLY:  average_2d, rotate_to_cosmo, find_horizontal_neighbours,&
                compute_horizontal_interp_weights,                              &
-               find_vertical_neighbours_and_weights, interpolate_2d,           &
-               gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, average_2d,  &
-               project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot
+               find_vertical_neighbours_and_weights_interp,                    &
+               find_vertical_neighbours_and_weights_average, interpolate_2d,   &
+               gamma_from_hemisphere, phic_to_phin, lamc_to_lamn,              &
+               project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot,   &
+               phirot2phi, rlarot2rla
     USE types
     USE util
@@ -80,4 +93,7 @@
     SAVE
     
+    REAL(dp) ::  averaging_angle   = 0.0_dp       !< latitudal and longitudal width of averaging regions [rad]
+    REAL(dp) ::  averaging_width_ns = 0.0_dp       !< longitudal width of averaging regions [m]
+    REAL(dp) ::  averaging_width_ew = 0.0_dp       !< latitudal width of averaging regions [m]
     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]
@@ -87,4 +103,10 @@
     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) ::  lam_centre        = 0.0_dp       !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  phi_centre        = 0.0_dp       !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  lam_east          = 0.0_dp       !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  lam_west          = 0.0_dp       !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  phi_north         = 0.0_dp       !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  phi_south         = 0.0_dp       !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [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]
@@ -100,8 +122,7 @@
     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) ::  f3                = 0.0_dp       !< Coriolis parameter
     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] 
@@ -140,7 +161,4 @@
     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 ::  average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging
-    INTEGER ::  k     !< indexing variable
     INTEGER ::  nt    !< number of output time steps
     INTEGER ::  nx    !< number of PALM-4U grid points in x direction
@@ -164,6 +182,4 @@
     LOGICAL ::  profile_grids_required
 
-    INTEGER ::  n_invar = 0 !< number of variables in the input variable table
-    INTEGER ::  n_outvar = 0 !< number of variables in the output variable table
     TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
     TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
@@ -219,8 +235,11 @@
     TYPE(grid_definition), TARGET ::  w_south_intermediate            !< 
     TYPE(grid_definition), TARGET ::  w_top_intermediate              !< 
-    TYPE(grid_definition), TARGET ::  scalar_profile_grid             !< 
-    TYPE(grid_definition), TARGET ::  scalar_profile_intermediate     !< 
-    TYPE(grid_definition), TARGET ::  w_profile_grid                  !< 
-    TYPE(grid_definition), TARGET ::  w_profile_intermediate          !< 
+
+    TYPE(grid_definition), TARGET ::  north_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  south_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  west_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  east_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  averaged_w_profile
 
     TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
@@ -237,4 +256,8 @@
     CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  radiation_files
     CHARACTER(LEN=SNAME) ::  input_prefix         !< Prefix of input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  flow_prefix          !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  soil_prefix          !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses
     CHARACTER(LEN=SNAME) ::  flow_suffix          !< Suffix of flow input files, e.g. 'flow'
     CHARACTER(LEN=SNAME) ::  soil_suffix          !< Suffix of soil input files, e.g. 'soil'
@@ -249,4 +272,5 @@
     
     SUBROUTINE setup_parameters()
+       INTEGER ::  k
 
 !
@@ -282,6 +306,6 @@
 
        ! Default atmospheric parameters
-       ug = 0.0_dp
-       vg = 0.0_dp
+       cfg % ug = 0.0_dp
+       cfg % vg = 0.0_dp
        cfg % p0 = P_SL
 
@@ -296,4 +320,5 @@
        input_prefix = 'laf'  ! 'laf' for COSMO-DE analyses
        cfg % flow_prefix = input_prefix
+       cfg % input_prefix = input_prefix
        cfg % soil_prefix = input_prefix
        cfg % radiation_prefix = input_prefix
@@ -304,4 +329,7 @@
        radiation_suffix = '-rad'
        soilmoisture_suffix = '-soilmoisture'
+
+       cfg % debug = .FALSE.
+       cfg % averaging_angle = 2.0_dp
 !
 !------------------------------------------------------------------------------
@@ -318,7 +346,17 @@
        cfg % ic_mode = 'volume'
        cfg % bc_mode = 'real'
-
-       ! Update default file names and domain centre
+       cfg % averaging_mode = 'level'
+
+       ! Overwrite defaults with user configuration
        CALL parse_command_line_arguments( cfg )
+
+       flow_prefix = TRIM(cfg % input_prefix)
+       radiation_prefix = TRIM(cfg % input_prefix)
+       soil_prefix = TRIM(cfg % input_prefix)
+       soilmoisture_prefix = TRIM(cfg % input_prefix)
+       IF (cfg % flow_prefix_is_set)  flow_prefix = TRIM(cfg % flow_prefix)
+       IF (cfg % radiation_prefix_is_set)  radiation_prefix = TRIM(cfg % radiation_prefix)
+       IF (cfg % soil_prefix_is_set)  soil_prefix = TRIM(cfg % soil_prefix)
+       IF (cfg % soilmoisture_prefix_is_set)  soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix)
 
        output_file % name = cfg % output_file
@@ -336,4 +374,6 @@
        END IF
 
+
+       ! Set default file paths, if not specified by user.
        CALL normalize_path(cfg % input_path)
        IF (TRIM(cfg % hhl_file) == '')  cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
@@ -361,12 +401,16 @@
 
        ! Generate input file lists
-       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)
+       CALL get_input_file_list(                                               &
+          cfg % start_date, start_hour_flow, end_hour, step_hour,              &
+          cfg % input_path, flow_prefix, flow_suffix, flow_files)
+       CALL get_input_file_list(                                               &
+          cfg % start_date, start_hour_soil, end_hour, step_hour,              &
+          cfg % input_path, soil_prefix, soil_suffix, soil_files)
+       CALL get_input_file_list(                                               &
+          cfg % start_date, start_hour_radiation, end_hour, step_hour,         &
+          cfg % input_path, 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, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
 
 !
@@ -499,4 +543,47 @@
 
  CALL run_control('time', 'comp')
+
+!------------------------------------------------------------------------------
+! Section 4.3: INIFOR averaging domains
+!------------------------------------------------------------------------------
+
+    ! Compute coordiantes of the PALM centre in the source (COSMO) system
+    phi_centre = phirot2phi(                                                   &
+       phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
+       rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
+       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
+       polgam = gam * TO_DEGREES                                               &
+    ) * TO_RADIANS
+
+    lam_centre = rlarot2rla(                                                   &
+       phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
+       rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
+       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
+       polgam = gam * TO_DEGREES                                               &
+    ) * TO_RADIANS
+
+    message = "PALM-4U centre:" // NEW_LINE('') // &
+       "           lon (lambda) = " // &
+       TRIM(real_to_str_f(lam_centre * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
+       "           lat (phi   ) = " // &
+       TRIM(real_to_str_f(phi_centre * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
+    CALL report( 'setup_parameters', message )
+
+
+    ! Compute boundaries of the central averaging box
+    averaging_angle = cfg % averaging_angle * TO_RADIANS
+    lam_east = lam_centre + 0.5_dp * averaging_angle
+    lam_west = lam_centre - 0.5_dp * averaging_angle
+    phi_north = phi_centre + 0.5_dp * averaging_angle
+    phi_south = phi_centre - 0.5_dp * averaging_angle
+    averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS
+    averaging_width_ns = averaging_angle * EARTH_RADIUS
+
+    ! Coriolis parameter
+    f3 = 2.0_dp * OMEGA * SIN(                                                 &
+       TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,&
+                              phi_n * TO_DEGREES, lambda_n * TO_DEGREES,       &
+                              gam * TO_DEGREES )                               &
+       )
 
     END SUBROUTINE setup_parameters
@@ -603,5 +690,5 @@
                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
                x0=x0, y0=y0, z0 = z0,                                          &
-               nx = nx, ny = ny, nz = nlev - 2)
+               nx = nx, ny = ny, nz = nlev - 1)
 
       IF (boundary_variables_required)  THEN
@@ -824,5 +911,5 @@
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2)
+                  nx = 0, ny = ny, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
@@ -830,5 +917,5 @@
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2)
+                  nx = 0, ny = ny, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
@@ -836,5 +923,5 @@
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2)
+                  nx = nx, ny = 0, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
@@ -842,5 +929,5 @@
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2)
+                  nx = nx, ny = 0, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
@@ -848,5 +935,5 @@
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny, nz = nlev - 2)
+                  nx = nx, ny = ny, nz = nlev - 1)
        END IF
 
@@ -856,32 +943,74 @@
 !------------------------------------------------------------------------------
 
-       profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR.              &
+       CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid,           &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='scalar')
+
+       CALL init_averaging_grid(averaged_w_profile, cosmo_grid,                &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0,              &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='w')
+
+       CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid,     &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_centre - averaging_angle,                          &
+               latmax = phi_centre,                                            &
+               kind='scalar')
+
+       CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid,     &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_centre,                                            &
+               latmax = phi_centre + averaging_angle,                          &
+               kind='scalar')
+
+       CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid,      &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_centre - averaging_angle,                          &
+               lonmax = lam_centre,                                            &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='scalar')
+
+       CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid,      &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_centre,                                            &
+               lonmax = lam_centre + averaging_angle,                          &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='scalar')
+
+       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,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  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,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  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,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  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,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nlev - 2, z=zw)
-       END IF
+       !IF (profile_grids_required)  THEN
+       !   ! Here, I use the 'boundary' kind, so the vertical coordinate uses the
+       !   ! PALM z coordinate.
+       !   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,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           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,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           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,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           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,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           nx = 0, ny = 0, nz = nlev - 2, z=zw)
+       !END IF
 
 !                                                                              
@@ -930,6 +1059,5 @@
 
        IF (TRIM(cfg % ic_mode) == 'profile')  THEN
-           CALL setup_averaging(cosmo_grid, palm_intermediate,                 &
-                average_imin, average_imax, average_jmin, average_jmax) 
+           !TODO: remove this conditional if not needed. 
        END IF
         
@@ -938,21 +1066,15 @@
 
 
-    SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode)
+    SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
 
        TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
-       TYPE(grid_definition), INTENT(INOUT), TARGET ::  palm_grid, palm_intermediate
+       TYPE(grid_definition), INTENT(INOUT), TARGET ::  grid, intermediate_grid
        CHARACTER, INTENT(IN)                        ::  kind
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
 
-       TYPE(grid_definition), POINTER      ::  grid
-       REAL(dp), DIMENSION(:), POINTER     ::  lat, lon
-       REAL(dp), DIMENSION(:,:,:), POINTER ::  h
-
-       LOGICAL :: setup_vertical
-
-       setup_vertical = .TRUE.
-       IF (PRESENT(ic_mode))  THEN
-          IF (TRIM(ic_mode) == 'profile')  setup_vertical = .FALSE.
-       END IF
+       REAL(dp), DIMENSION(:), POINTER     ::  cosmo_lat, cosmo_lon
+       REAL(dp), DIMENSION(:,:,:), POINTER ::  cosmo_h
+
+       LOGICAL :: setup_volumetric
 
 !------------------------------------------------------------------------------
@@ -964,25 +1086,25 @@
        CASE('s') ! scalars
 
-          lat => cosmo_grid % lat
-          lon => cosmo_grid % lon
-          h   => cosmo_grid % hfl
+          cosmo_lat => cosmo_grid % lat
+          cosmo_lon => cosmo_grid % lon
+          cosmo_h   => cosmo_grid % hfl
 
        CASE('w') ! vertical velocity
 
-          lat => cosmo_grid % lat
-          lon => cosmo_grid % lon
-          h   => cosmo_grid % hhl
+          cosmo_lat => cosmo_grid % lat
+          cosmo_lon => cosmo_grid % lon
+          cosmo_h   => cosmo_grid % hhl
 
        CASE('u') ! x velocity
 
-          lat => cosmo_grid % lat
-          lon => cosmo_grid % lonu
-          h   => cosmo_grid % hfl
+          cosmo_lat => cosmo_grid % lat
+          cosmo_lon => cosmo_grid % lonu
+          cosmo_h   => cosmo_grid % hfl
 
        CASE('v') ! y velocity
 
-          lat => cosmo_grid % latv
-          lon => cosmo_grid % lon
-          h   => cosmo_grid % hfl
+          cosmo_lat => cosmo_grid % latv
+          cosmo_lon => cosmo_grid % lon
+          cosmo_h   => cosmo_grid % hfl
 
        CASE DEFAULT
@@ -993,11 +1115,12 @@
        END SELECT
 
-       grid => palm_intermediate
-
-       CALL find_horizontal_neighbours(lat, lon,                               &
-          grid % clat, grid % clon, grid % ii, grid % jj)
-
-       CALL compute_horizontal_interp_weights(lat, lon,                        &
-          grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz)
+       CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon,                   &
+          intermediate_grid % clat, intermediate_grid % clon,                  &
+          intermediate_grid % ii, intermediate_grid % jj)
+
+       CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,            &
+          intermediate_grid % clat, intermediate_grid % clon,                  &
+          intermediate_grid % ii, intermediate_grid % jj,                      &
+          intermediate_grid % w_horiz)
 
 !------------------------------------------------------------------------------
@@ -1005,55 +1128,27 @@
 !------------------------------------------------------------------------------
 
-       IF (setup_vertical)  THEN
-          ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 
-          grid % h(:,:,:) = - EARTH_RADIUS
+       ! If profile initialization is chosen, we--somewhat counterintuitively--
+       ! don't need to compute vertical interpolation weights. At least, we
+       ! don't need them on the intermediate grid, which fills the entire PALM
+       ! domain volume. Instead we need vertical weights for the intermediate
+       ! profile grids, which get computed in setup_averaging().
+       setup_volumetric = .TRUE.
+       IF (PRESENT(ic_mode))  THEN
+          IF (TRIM(ic_mode) == 'profile')  setup_volumetric = .FALSE.
+       END IF
+
+       IF (setup_volumetric)  THEN
+          ALLOCATE( intermediate_grid % h(0:intermediate_grid % nx,            &
+                                          0:intermediate_grid % ny,            &
+                                          0:intermediate_grid % nz) ) 
+          intermediate_grid % h(:,:,:) = - EARTH_RADIUS
 
           ! For w points, use hhl, for scalars use hfl
           ! compute the full heights for the intermediate grids
-          CALL interpolate_2d(cosmo_grid % hfl, grid % h, grid)
-          CALL find_vertical_neighbours_and_weights(palm_grid, grid)
+          CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid)
+          CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
        END IF
        
     END SUBROUTINE setup_interpolation
-
-
-    SUBROUTINE setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax)
-
-       TYPE(grid_definition), INTENT(IN) ::  cosmo_grid, palm_intermediate
-       INTEGER, INTENT(INOUT)            ::  imin, imax, jmin, jmax
-
-       TYPE(grid_definition), POINTER    ::  grid
-       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_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)
-       imax = CEILING(lonmax_pos)
-       jmin = FLOOR(latmin_pos)
-       jmax = CEILING(latmax_pos)
-       
-       ! average heights for intermediate scalar and w profile grids
-       grid => scalar_profile_intermediate
-       ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 
-       grid % h(:,:,:) = - EARTH_RADIUS
-       CALL average_2d(cosmo_grid % hfl, grid % h(0,0,:), imin, imax, jmin, jmax)
-       CALL find_vertical_neighbours_and_weights(scalar_profile_grid, grid)
-
-       grid => w_profile_intermediate
-       ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 
-       grid % h(:,:,:) = - EARTH_RADIUS
-       CALL average_2d(cosmo_grid % hhl, grid % h(0,0,:), imin, imax, jmin, jmax)
-       CALL find_vertical_neighbours_and_weights(w_profile_grid, grid)
-
-    END SUBROUTINE setup_averaging
-
 
 !------------------------------------------------------------------------------!
@@ -1291,4 +1386,177 @@
 
     END SUBROUTINE init_grid_definition
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initializes averagin_grid-type variables.
+!> 
+!> Input parameters:
+!> -----------------
+!>
+!> cosmo_grid : grid_definition-type varieble initialized as 'cosmo-de' grid
+!>    providing longitudes and latitudes of the parent grid
+!>
+!> x, y : location of the profile to be averaged in the PALM system [m]
+!>
+!> z : vertical grid coordinates of the profile in the PALM system [m]
+!>
+!> <lat/lon><min/max>: boundaries of the averaging region in the parent system,
+!>    i.e. the COSMO-DE rotated-pole system. [rad]
+!>
+!> kind : kind of quantity to be averaged using this averaging grid.
+!>    Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the
+!>    parent/COSMO columns for the region in get_cosmo_averaging_region() is
+!>    independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity
+!>    components, which are computed in the preprocessing step.
+!>
+!> Output parameters:
+!> ------------------
+!> avg_grid : averaging_grid variable to be initialized.
+!------------------------------------------------------------------------------!
+    SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0,          &
+       lonmin, lonmax, latmin, latmax, kind)
+
+       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
+       TYPE(grid_definition), INTENT(IN)    ::  cosmo_grid
+       REAL(dp), INTENT(IN)                 ::  x, y, z0
+       REAL(dp), INTENT(IN), TARGET         ::  z(:)
+       REAL(dp), INTENT(IN)                 ::  lonmin, lonmax, latmin, latmax
+
+       CHARACTER(LEN=*), INTENT(IN)         ::  kind
+
+       LOGICAL                              :: level_based_averaging
+
+       INTEGER :: i, j
+
+       ALLOCATE( avg_grid % x(1) )
+       ALLOCATE( avg_grid % y(1) )
+       avg_grid % x(1) = x
+       avg_grid % y(1) = y
+       avg_grid % z => z
+       avg_grid % z0 = z0
+
+       avg_grid % nz = SIZE(z, 1)
+
+       ALLOCATE( avg_grid % lon(2) )
+       ALLOCATE( avg_grid % lat(2) )
+       avg_grid % lon(1:2) = (/lonmin, lonmax/)
+       avg_grid % lat(1:2) = (/latmin, latmax/)
+
+       avg_grid % kind = TRIM(kind)
+
+       ! Find and store COSMO columns that fall into the coordinate range
+       ! given by avg_grid % clon, % clat
+       CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
+
+       ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) )
+       ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) )
+       ! Compute average COSMO levels in the averaging region
+       SELECT CASE(avg_grid % kind)
+
+       CASE('scalar', 'u', 'v')
+          avg_grid % cosmo_h => cosmo_grid % hfl
+
+       CASE('w')
+          avg_grid % cosmo_h => cosmo_grid % hhl
+
+       CASE DEFAULT
+          message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
+                    "' is not supported. Use 'scalar', 'u', or 'v'."
+          CALL abort('get_cosmo_averaging_region', message)
+
+       END SELECT
+
+       ! For level-besed averaging, compute average heights
+       level_based_averaging = .TRUE.
+       IF (level_based_averaging)  THEN
+          ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) )
+ 
+          CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:),                            &
+                          avg_grid % iii, avg_grid % jjj)
+       END IF
+
+       ! Copy averaged grid to all COSMO columnts, leads to computing the same
+       ! vertical interpolation weights for all columns and to COSMO grid level
+       ! based averaging onto the averaged COSMO heights.
+       IF ( TRIM(cfg % averaging_mode) == 'level' )  THEN
+          FORALL(                                                              &
+             i=1:SIZE(avg_grid % cosmo_h, 1),                                  &
+             j=1:SIZE(avg_grid % cosmo_h, 2)                                   &
+          )   avg_grid % cosmo_h(i,j,:) = avg_grid % h(1,1,:)
+       END IF
+
+       ! Compute vertical weights and neighbours
+       CALL find_vertical_neighbours_and_weights_average(avg_grid)
+
+    END SUBROUTINE init_averaging_grid
+
+
+    SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
+       TYPE(grid_definition), INTENT(INOUT)         ::  avg_grid
+       TYPE(grid_definition), TARGET, INTENT(IN)    ::  cosmo_grid
+
+       REAL(dp), DIMENSION(:), POINTER              ::  cosmo_lon, cosmo_lat
+       REAL(dp)                                     ::  dlon, dlat
+
+       INTEGER ::  i, j, imin, imax, jmin, jmax, l, nx, ny
+
+
+       SELECT CASE( TRIM(avg_grid % kind) )
+
+       CASE('scalar', 'w')
+          cosmo_lon => cosmo_grid % lon
+          cosmo_lat => cosmo_grid % lat
+
+       CASE('u')
+          cosmo_lon => cosmo_grid % lonu
+          cosmo_lat => cosmo_grid % lat
+
+       CASE('v')
+          cosmo_lon => cosmo_grid % lon
+          cosmo_lat => cosmo_grid % latv
+
+       CASE DEFAULT
+          message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
+                    "' is not supported. Use 'scalar', 'u', or 'v'."
+          CALL abort('get_cosmo_averaging_region', message)
+
+       END SELECT
+
+
+       ! FIXME: make dlon, dlat parameters of the grid_defintion type
+       dlon = cosmo_lon(1) - cosmo_lon(0)
+       dlat = cosmo_lat(1) - cosmo_lat(0)
+
+       imin = CEILING( (avg_grid % lon(1) - cosmo_lon(0)) / dlon )
+       imax = FLOOR(   (avg_grid % lon(2) - cosmo_lon(0)) / dlon )
+
+       jmin = CEILING( (avg_grid % lat(1) - cosmo_lat(0)) / dlat )
+       jmax = FLOOR(   (avg_grid % lat(2) - cosmo_lat(0)) / dlat )
+       
+       message = "Averaging over "//                                           &
+                 TRIM(str(imin)) // " < i < " // TRIM(str(imax)) //            &
+                 " and " //                                                    &
+                 TRIM(str(jmin)) // " < j < " // TRIM(str(jmax))
+       CALL report( 'get_cosmo_averaging_region', message )
+
+       nx = imax - imin + 1
+       ny = jmax - jmin + 1
+       avg_grid % n_columns = nx * ny
+
+       ALLOCATE( avg_grid % iii(avg_grid % n_columns),                                 &
+                 avg_grid % jjj(avg_grid % n_columns) )
+
+       l = 0
+       DO j = jmin, jmax
+       DO i = imin, imax
+          l = l + 1
+          avg_grid % iii(l) = i
+          avg_grid % jjj(l) = j
+       END DO
+       END DO
+
+    END SUBROUTINE get_cosmo_averaging_region
 
 
@@ -1718,25 +1986,30 @@
        )
 
-       !potential temperature, surface pressure, including nudging and subsidence
+       !potential temperature, surface pressure, specific humidity including
+       !nudging and subsidence, and geostrophic winds ug, vg
        io_group_list(3) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          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 including nudging and subsidence
-       io_group_list(4) = init_io_group(                                       &
-          in_files = flow_files,                                               &
-          out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
-          in_var_list = input_var_table(4:4),                                  &
-          kind = 'scalar'                                                      &
-       )
-
-       !u and v velocity and geostrophic winds ug, vg
+          out_vars = [output_var_table(56:64),                                 & ! internal averaged density and pressure profiles
+                      output_var_table(3:8), output_var_table(9:14),           &
+                      output_var_table(42:42), output_var_table(43:44),        &
+                      output_var_table(49:51), output_var_table(52:54)],       &
+          in_var_list = (/input_var_table(3), input_var_table(17),             & ! T, P, QV
+                          input_var_table(4) /),                               &
+          kind = 'thermodynamics',                                             &
+          n_output_quantities = 4                                              & ! P, Theta, Rho, qv
+       )
+
+       !Moved to therodynamic io_group
+       !io_group_list(4) = init_io_group(                                       &
+       !   in_files = flow_files,                                               &
+       !   out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
+       !   in_var_list = input_var_table(4:4),                                  &
+       !   kind = 'scalar'                                                      &
+       !)
+
+       !u and v velocity, ug anv vg moved to thermodynamic io group. 
        io_group_list(5) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          out_vars = [output_var_table(15:26), output_var_table(43:46)],       &
+          out_vars = [output_var_table(15:26), output_var_table(45:46)],       &
           !out_vars = output_var_table(15:20),                                  &
           in_var_list = input_var_table(5:6),                                  &
@@ -1845,9 +2118,11 @@
 
 
-    FUNCTION init_io_group(in_files, out_vars, in_var_list, kind) RESULT(group)
+    FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
+                           n_output_quantities) RESULT(group)
        CHARACTER(LEN=PATH), INTENT(IN) ::  in_files(:)
        CHARACTER(LEN=*), INTENT(IN)    ::  kind
        TYPE(nc_var), INTENT(IN)        ::  out_vars(:)
        TYPE(nc_var), INTENT(IN)        ::  in_var_list(:)
+       INTEGER, OPTIONAL               ::  n_output_quantities
 
        TYPE(io_group)                  ::  group
@@ -1855,7 +2130,14 @@
        group % nt = SIZE(in_files)
        group % nv = SIZE(out_vars)
+       group % n_inputs = SIZE(in_var_list)
        group % kind = TRIM(kind)
 
-       ALLOCATE(group % in_var_list( SIZE(in_var_list) ))
+       IF ( PRESENT(n_output_quantities) )  THEN
+          group % n_output_quantities = n_output_quantities
+       ELSE
+          group % n_output_quantities = group % n_inputs
+       END IF
+
+       ALLOCATE(group % in_var_list(group % n_inputs))
        ALLOCATE(group % in_files(group % nt))
        ALLOCATE(group % out_vars(group % nv))
@@ -1876,5 +2158,5 @@
     SUBROUTINE fini_grids()
 
-       CALL report('fini_grids', 'Deallocating grids')
+       CALL report('fini_grids', 'Deallocating grids', cfg % debug)
        
        DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
@@ -1904,4 +2186,6 @@
     SUBROUTINE setup_variable_tables(ic_mode)
        CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
+       INTEGER                      ::  n_invar = 0  !< number of variables in the input variable table
+       INTEGER                      ::  n_outvar = 0 !< number of variables in the output variable table
        TYPE(nc_var), POINTER        ::  var
 
@@ -1914,5 +2198,5 @@
 
        n_invar = 17
-       n_outvar = 55
+       n_outvar = 64
        ALLOCATE( input_var_table(n_invar) )
        ALLOCATE( output_var_table(n_outvar) )
@@ -2052,6 +2336,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(3) % grid => scalar_profile_grid
-          output_var_table(3) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(3) % averaging_grid => averaged_scalar_profile
        END IF
 
@@ -2122,5 +2405,5 @@
           units             = "kg/kg",                                         &
           kind              = "init scalar",                                   &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = palm_grid,                                       &
@@ -2129,6 +2412,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(9) % grid => scalar_profile_grid
-          output_var_table(9) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(9) % averaging_grid => averaged_scalar_profile
        END IF
 
@@ -2139,5 +2421,5 @@
           units             = "kg/kg",                                         &
           kind              = "left scalar",                                   &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_west_grid,                               &
@@ -2151,5 +2433,5 @@
           units             = "kg/kg",                                         &
           kind              = "right scalar",                                  &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_east_grid,                               &
@@ -2163,5 +2445,5 @@
           units             = "kg/kg",                                         &
           kind              = "north scalar",                                  &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_north_grid,                              &
@@ -2175,5 +2457,5 @@
           units             = "kg/kg",                                         &
           kind              = "south scalar",                                  &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_south_grid,                              &
@@ -2187,5 +2469,5 @@
           units             = "kg/kg",                                         &
           kind              = "top scalar",                                    &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_top_grid,                                &
@@ -2206,6 +2488,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(15) % grid => scalar_profile_grid
-          output_var_table(15) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(15) % averaging_grid => averaged_scalar_profile
        END IF 
 
@@ -2283,6 +2564,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(21) % grid => scalar_profile_grid
-          output_var_table(21) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(21) % averaging_grid => averaged_scalar_profile
        END IF
 
@@ -2360,6 +2640,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(27) % grid => w_profile_grid
-          output_var_table(27) % intermediate_grid => w_profile_intermediate
+          output_var_table(27) % averaging_grid => averaged_w_profile
        END IF
 
@@ -2546,4 +2825,5 @@
           intermediate_grid = palm_intermediate                                &
        )
+       output_var_table(42) % averaging_grid => averaged_scalar_profile
 
        output_var_table(43) = init_nc_var(                                     &
@@ -2552,9 +2832,9 @@
           long_name         = "geostrophic wind (u component)",                &
           units             = "m/s",                                           &
-          kind              = "constant scalar profile",                       &
+          kind              = "geostrophic",                                   &
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
 
@@ -2564,9 +2844,9 @@
           long_name         = "geostrophic wind (v component)",                &
           units             = "m/s",                                           &
-          kind              = "constant scalar profile",                       &
+          kind              = "geostrophic",                                   &
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
 
@@ -2576,9 +2856,9 @@
           long_name         = "wind component in x direction",                 &
           units             = "m/s",                                           &
-          kind              = "large-scale scalar forcing",                    &
+          kind              = "geostrophic",                                   &
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(45) % to_be_processed = ls_forcing_variables_required
@@ -2592,6 +2872,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(46) % to_be_processed = ls_forcing_variables_required
@@ -2605,6 +2885,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = w_profile_grid,                                  &
-          intermediate_grid = w_profile_intermediate                           &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(47) % to_be_processed = ls_forcing_variables_required
@@ -2618,6 +2898,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = w_profile_grid,                                  &
-          intermediate_grid = w_profile_intermediate                           &
+          grid              = averaged_w_profile,                              &
+          intermediate_grid = averaged_w_profile                               &
        )
        output_var_table(48) % to_be_processed = ls_forcing_variables_required
@@ -2632,6 +2912,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(49) % to_be_processed = ls_forcing_variables_required
@@ -2645,6 +2925,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(50) % to_be_processed = ls_forcing_variables_required
@@ -2658,6 +2938,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(51) % to_be_processed = ls_forcing_variables_required
@@ -2669,8 +2949,8 @@
           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                      &
+          input_id          = 3,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(52) % to_be_processed = ls_forcing_variables_required
@@ -2683,8 +2963,8 @@
           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                      &
+          input_id          = 3,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(53) % to_be_processed = ls_forcing_variables_required
@@ -2696,8 +2976,8 @@
           units             = "kg/kg",                                         &
           kind              = "large-scale scalar forcing",                    &
-          input_id          = 1,                                               &
-          output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          input_id          = 3,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(54) % to_be_processed = ls_forcing_variables_required
@@ -2711,9 +2991,141 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(55) % to_be_processed = ls_forcing_variables_required
 
+
+       output_var_table(56) = init_nc_var(                                     &
+          name              = 'internal_density_centre',                              &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
+       )
+       output_var_table(56) % averaging_grid => averaged_scalar_profile
+
+
+       output_var_table(57) = init_nc_var(                                     &
+          name              = 'internal_density_north',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = north_averaged_scalar_profile,                   &
+          intermediate_grid = north_averaged_scalar_profile                    &
+       )
+       output_var_table(57) % averaging_grid => north_averaged_scalar_profile
+       output_var_table(57) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(58) = init_nc_var(                                     &
+          name              = 'internal_density_south',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = south_averaged_scalar_profile,                   &
+          intermediate_grid = south_averaged_scalar_profile                    &
+       )
+       output_var_table(58) % averaging_grid => south_averaged_scalar_profile
+       output_var_table(58) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(59) = init_nc_var(                                     &
+          name              = 'internal_density_east',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = east_averaged_scalar_profile,                    &
+          intermediate_grid = east_averaged_scalar_profile                     &
+       )
+       output_var_table(59) % averaging_grid => east_averaged_scalar_profile
+       output_var_table(59) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(60) = init_nc_var(                                     &
+          name              = 'internal_density_west',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = west_averaged_scalar_profile,                    &
+          intermediate_grid = west_averaged_scalar_profile                     &
+       )
+       output_var_table(60) % averaging_grid => west_averaged_scalar_profile
+       output_var_table(60) % to_be_processed = .NOT. cfg % ug_is_set
+
+       output_var_table(61) = init_nc_var(                                     &
+          name              = 'internal_pressure_north',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = north_averaged_scalar_profile,                   &
+          intermediate_grid = north_averaged_scalar_profile                    &
+       )
+       output_var_table(61) % averaging_grid => north_averaged_scalar_profile
+       output_var_table(61) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(62) = init_nc_var(                                     &
+          name              = 'internal_pressure_south',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = south_averaged_scalar_profile,                   &
+          intermediate_grid = south_averaged_scalar_profile                    &
+       )
+       output_var_table(62) % averaging_grid => south_averaged_scalar_profile
+       output_var_table(62) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(63) = init_nc_var(                                     &
+          name              = 'internal_pressure_east',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = east_averaged_scalar_profile,                    &
+          intermediate_grid = east_averaged_scalar_profile                     &
+       )
+       output_var_table(63) % averaging_grid => east_averaged_scalar_profile
+       output_var_table(63) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(64) = init_nc_var(                                     &
+          name              = 'internal_pressure_west',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = west_averaged_scalar_profile,                    &
+          intermediate_grid = west_averaged_scalar_profile                     &
+       )
+       output_var_table(64) % averaging_grid => west_averaged_scalar_profile
+       output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set
 
        ! Attributes shared among all variables
@@ -2732,6 +3144,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
@@ -2772,4 +3184,5 @@
           var % dimvarids(1:3)  = output_file % dimvarids_soil
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_2d"
 
@@ -2781,4 +3194,5 @@
           var % dimvarids(1:3)  = output_file % dimvarids_scl
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2794,4 +3208,5 @@
           var % dimvarids(3)    = output_file % dimvarids_scl(3)
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2807,4 +3222,5 @@
           var % dimvarids(3)    = output_file % dimvarids_scl(3)
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2820,4 +3236,5 @@
           var % dimvarids(3)    = output_file % dimvarids_vel(3)
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2829,4 +3246,5 @@
           var % dimvarids(1)    = output_file % dimvarids_scl(3) !z
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average profile"
 
@@ -2838,4 +3256,5 @@
           var % dimvarids(1)    = output_file % dimvarids_vel(3) !z
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average profile"
 
@@ -2848,4 +3267,5 @@
           var % dimvarids(1:2)  = output_file % dimvarids_soil(1:2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_2d"
 
@@ -2860,4 +3280,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2872,4 +3293,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2884,4 +3306,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2896,4 +3319,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2908,4 +3332,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2920,4 +3345,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2932,4 +3358,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2944,4 +3371,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2956,4 +3384,5 @@
           var % dimvarids(2)    = output_file % dimvarids_vel(2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2968,4 +3397,5 @@
           var % dimvarids(2)    = output_file % dimvarids_vel(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2980,4 +3410,5 @@
           var % dimvarids(2)    = output_file % dimvarids_vel(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2988,5 +3419,6 @@
           var % dimvarids(1)    = output_file % dimvarid_time
           var % to_be_processed = .TRUE.
-          var % task            = "average scalar"
+          var % is_internal     = .FALSE.
+          var % task            = "average profile"
 
        CASE( 'constant scalar profile' )
@@ -2998,4 +3430,5 @@
           var % dimvarids(1)    = output_file % dimvarids_scl(3)
           var % to_be_processed = .TRUE.
+          var % is_internal     = .FALSE.
           var % task            = "set profile"
 
@@ -3008,5 +3441,17 @@
           var % dimvarids(1)    = output_file % dimvarids_scl(3)
           var % to_be_processed = ls_forcing_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average large-scale profile"
+
+       CASE( 'geostrophic' )
+          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 = .TRUE.
+          var % is_internal     = .FALSE.
+          var % task            = "geostrophic winds"
 
        CASE( 'large-scale w forcing' )
@@ -3018,5 +3463,17 @@
           var % dimvarids(1)    = output_file % dimvarids_vel(3)
           var % to_be_processed = ls_forcing_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average large-scale profile"
+
+       CASE( 'internal profile' )
+          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 = .TRUE.
+          var % is_internal     = .TRUE.
+          var % task            = "internal profile"
 
        CASE DEFAULT
@@ -3031,5 +3488,5 @@
     SUBROUTINE fini_variables()
 
-       CALL report('fini_variables', 'Deallocating variable table')
+       CALL report('fini_variables', 'Deallocating variable table', cfg % debug)
        DEALLOCATE( input_var_table )
 
@@ -3039,5 +3496,5 @@
     SUBROUTINE fini_io_groups()
 
-       CALL report('fini_io_groups', 'Deallocating IO groups')
+       CALL report('fini_io_groups', 'Deallocating IO groups', cfg % debug)
        DEALLOCATE( io_group_list )
 
@@ -3047,5 +3504,5 @@
     SUBROUTINE fini_file_lists()
        
-       CALL report('fini_file_lists', 'Deallocating file lists')
+       CALL report('fini_file_lists', 'Deallocating file lists', cfg % debug)
        DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files )
 
@@ -3134,6 +3591,6 @@
              ! rotated velocities
              DO k = 1, nz
-             DO j = 2, ny
-             DO i = 2, nx
+             DO j = 1, ny
+             DO i = 1, nx
                 CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k),     &
                                vrot = preprocess_buffer(2) % array(i,j,k),     &
@@ -3171,5 +3628,5 @@
           CALL report('preprocess', message)
        
-       CASE( 'temperature' ) ! P and T
+       CASE( 'thermodynamics' ) ! T, P, QV
           nx = SIZE(input_buffer(1) % array, 1)
           ny = SIZE(input_buffer(1) % array, 2)
@@ -3192,5 +3649,6 @@
 
                 input_buffer (2) % array(i,j,:) =                              &
-                   input_buffer (2) % array(i,j,:) + basic_state_pressure(:)
+                   HECTO * input_buffer (2) % array(i,j,:) +                   &
+                   basic_state_pressure(:)
 
              END DO
@@ -3204,10 +3662,42 @@
 
           END IF
+          ! pressure
+          input_buffer(2) % is_preprocessed = .TRUE.
+
+          ! Copy temperature to last input buffer array
+          ALLOCATE(                                                            &
+              input_buffer( group % n_output_quantities ) % array (nx, ny, nz) &
+          )
+          input_buffer(group % n_output_quantities) % array(:,:,:) =           &
+              input_buffer(1) % array(:,:,:)
 
           ! Convert absolute to potential temperature
-          input_buffer(1) % array(:,:,:) = input_buffer(1) % array(:,:,:) *    &
-             EXP( RD_PALM/CP_PALM * LOG( P_REF / input_buffer(2) % array(:,:,:) ))
-
-          input_buffer(1:2) % is_preprocessed = .TRUE.
+          CALL potential_temperature(                                          &
+             t = input_buffer(1) % array(:,:,:),                               &
+             p = input_buffer(2) % array(:,:,:),                               &
+             p_ref = P_REF,                                                    &
+             r = RD_PALM,                                                      &
+             cp = CP_PALM                                                      &
+          )
+
+          ! potential temperature
+          input_buffer(1) % is_preprocessed = .TRUE.
+
+          ! Convert temperature copy to density
+          CALL moist_density(                                                  &
+             t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), &
+             p = input_buffer(2) % array(:,:,:),                               &
+             qv = input_buffer(3) % array(:,:,:),                              &
+             rd = RD,                                                          &
+             rv = RV                                                           &
+          )
+
+          ! qv
+          input_buffer(3) % is_preprocessed = .TRUE.
+
+          ! density
+          input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE.
+
+
           message = "Input buffers for group '" // TRIM(group % kind) // "'"//&
              " preprocessed sucessfully."
@@ -3336,5 +3826,5 @@
 
        CASE DEFAULT
-          message = "Buffer kind '" // TRIM(group % kind) // "' is not supported."
+          message = "IO group kind '" // TRIM(group % kind) // "' is not supported."
           CALL abort('prerpocess', message)
 
Index: /palm/trunk/UTIL/inifor/src/inifor.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/inifor.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/inifor.f90	(revision 3395)
@@ -26,4 +26,9 @@
 ! -----------------
 ! $Id$
+! Added main loop support for computation of geostrophic winds and surface
+!     pressure
+! Cleaned up terminal output, show some meVssages only with --debug option
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Introduced new PALM grid stretching
 ! Renamend initial-condition mode variable 'mode' to 'ic_mode'
@@ -56,10 +61,12 @@
                fini_file_lists, preprocess, origin_lon, origin_lat,            &
                output_file, io_group_list, output_var_table,                   &
-               cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, cfg,             &
-               average_imin, average_imax, average_jmin, average_jmax
-
+               cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3,                 &
+               averaging_width_ns, averaging_width_ew, phi_n, lambda_n,        &
+               lam_centre, phi_centre
     USE io
     USE transform,                                                             &
-        ONLY:  average_profile, average_2d, interpolate_2d, interpolate_3d
+        ONLY:  average_profile, interpolate_2d, interpolate_3d,                &
+               geostrophic_winds, extrapolate_density, extrapolate_pressure,   &
+               get_surface_pressure
     USE types
     
@@ -70,7 +77,22 @@
     INTEGER                                 ::  iter
     REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) ::  output_arr
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_arr
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_arr
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west
+    REAL(dp), POINTER, DIMENSION(:)         ::  internal_arr
+    REAL(dp), POINTER, DIMENSION(:)         ::  ug_vg_arr
     TYPE(nc_var), POINTER                   ::  output_var
     TYPE(io_group), POINTER                 ::  group
     TYPE(container), ALLOCATABLE            ::  input_buffer(:)
+    LOGICAL, SAVE                           ::  ug_vg_have_been_computed = .FALSE.
+    LOGICAL, SAVE                           ::  debugging_output = .TRUE.
     
 !> \mainpage About INIFOR
@@ -199,12 +221,152 @@
  CALL run_control('time', 'alloc')
                       
-
                       CALL average_profile(                                    &
                          input_buffer(output_var % input_id) % array(:,:,:),   &
-                         output_arr(:,:,:), average_imin, average_imax,        &
-                         average_jmin, average_jmax,                           &
-                         output_var % intermediate_grid,                       &
-                         output_var % grid)
- CALL run_control('time', 'comp')
+                         output_arr(0,0,:),                                    &
+                         output_var % averaging_grid)
+
+                      IF ( TRIM(output_var % name) ==                          &
+                           'surface_forcing_surface_pressure' )  THEN
+
+                         IF ( cfg % p0_is_set )  THEN
+                            output_arr(0,0,1) = p0
+                         ELSE
+                            CALL get_surface_pressure(                         &
+                               output_arr(0,0,:), rho_centre,                  &
+                               output_var % averaging_grid)
+                         END IF
+
+                      END IF
+ CALL run_control('time', 'comp')
+
+                   CASE ( 'internal profile' )
+
+                      message = "Averaging of internal profile for variable '" //&
+                         TRIM(output_var % name) // "' is not supported."
+
+                      SELECT CASE (TRIM(output_var % name))
+
+                      CASE('internal_density_centre')
+                         ALLOCATE( rho_centre( 1:output_var % grid % nz ) )
+                         internal_arr => rho_centre
+
+                      CASE('internal_density_north')
+                         ALLOCATE( rho_north( 1:output_var % grid % nz ) )
+                         internal_arr => rho_north
+
+                      CASE('internal_density_south')
+                         ALLOCATE( rho_south( 1:output_var % grid % nz ) )
+                         internal_arr => rho_south
+
+                      CASE('internal_density_east')
+                         ALLOCATE( rho_east( 1:output_var % grid % nz) )
+                         internal_arr => rho_east
+
+                      CASE('internal_density_west')
+                         ALLOCATE( rho_west( 1:output_var % grid % nz ) )
+                         internal_arr => rho_west
+
+                      CASE('internal_pressure_north')
+                         ALLOCATE( p_north( 1:output_var % grid % nz ) )
+                         internal_arr => p_north
+
+                      CASE('internal_pressure_south')
+                         ALLOCATE( p_south( 1:output_var % grid % nz ) )
+                         internal_arr => p_south
+
+                      CASE('internal_pressure_east')
+                         ALLOCATE( p_east( 1:output_var % grid % nz) )
+                         internal_arr => p_east
+
+                      CASE('internal_pressure_west')
+                         ALLOCATE( p_west( 1:output_var % grid % nz ) )
+                         internal_arr => p_west
+
+                      CASE DEFAULT
+                         CALL abort('main loop', message)
+
+                      END SELECT
+ CALL run_control('time', 'alloc')
+
+
+                      CALL average_profile(                                 &
+                         input_buffer(output_var % input_id) % array(:,:,:),&
+                         internal_arr(:),                                   &
+                         output_var % averaging_grid)
+
+                      SELECT CASE (TRIM(output_var % name))
+
+                      CASE('internal_density_centre',                          &
+                           'internal_density_north',                           &
+                           'internal_density_south',                           &
+                           'internal_density_east',                            &
+                           'internal_density_west')
+                         CALL extrapolate_density(internal_arr,                &
+                                                  output_var % averaging_grid)
+
+                      CASE('internal_pressure_north')
+                         CALL extrapolate_pressure(internal_arr, rho_north,    &
+                                                   output_var % averaging_grid)
+
+                      CASE('internal_pressure_south')
+                         CALL extrapolate_pressure(internal_arr, rho_south,    &
+                                                   output_var % averaging_grid)
+
+                      CASE('internal_pressure_east')
+                         CALL extrapolate_pressure(internal_arr, rho_east,     &
+                                                   output_var % averaging_grid)
+
+                      CASE('internal_pressure_west')
+                         CALL extrapolate_pressure(internal_arr, rho_west,     &
+                                                   output_var % averaging_grid)
+
+                      CASE DEFAULT
+                         CALL abort('main loop', message)
+
+                      END SELECT
+
+                      IF (.TRUE.)  THEN
+                         ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
+                         output_arr(1,1,:) = internal_arr(:)
+                      END IF
+ CALL run_control('time', 'comp')
+
+
+                   ! This case gets called twice, the first time for ug, the
+                   ! second time for vg. We compute ug and vg at the first call
+                   ! and keep vg (and ug for that matter) around for the second
+                   ! call.
+                   CASE ( 'geostrophic winds' )
+
+                      IF (.NOT. ug_vg_have_been_computed )  THEN
+                         ALLOCATE( ug_arr(output_var % grid % nz) )
+                         ALLOCATE( vg_arr(output_var % grid % nz) )
+
+                         IF ( cfg % ug_is_set )  THEN
+                            ug_arr = cfg % ug
+                            vg_arr = cfg % vg
+                         ELSE
+                            CALL geostrophic_winds( p_north, p_south, p_east,     &
+                                                    p_west, rho_centre, f3,      &
+                                                    averaging_width_ew,           &
+                                                    averaging_width_ns,           &
+                                                    phi_n, lambda_n, phi_centre, lam_centre, &
+                                                    ug_arr, vg_arr )
+                         END IF
+
+                         ug_vg_have_been_computed = .TRUE.
+
+                      END IF
+
+                      ! Prepare output of geostrophic winds
+                      SELECT CASE(TRIM(output_var % name))
+                      CASE ('ls_forcing_ug')
+                         ug_vg_arr => ug_arr
+                      CASE ('ls_forcing_vg')
+                         ug_vg_arr => vg_arr
+                      END SELECT
+
+                      ALLOCATE( output_arr(1,1,output_var % grid % nz) )
+                      output_arr(1,1,:) = ug_vg_arr(:)
 
                    CASE ( 'average scalar' )
@@ -222,9 +384,9 @@
                       SELECT CASE (TRIM(output_var % name))
 
-                      CASE('ls_forcing_ug')
-                          output_arr(1, 1, :) = ug
-
-                      CASE('ls_forcing_vg')
-                          output_arr(1, 1, :) = vg
+                      !CASE('ls_forcing_ug')
+                      !    output_arr(1, 1, :) = ug
+
+                      !CASE('ls_forcing_vg')
+                      !    output_arr(1, 1, :) = vg
 
                       CASE('nudging_tau')
@@ -256,10 +418,12 @@
 !- Section 2.3: Write current time step of current variable
 !------------------------------------------------------------------------------
-                   message = "Writing variable '" // TRIM(output_var%name) // "'."
-                   CALL report('main loop', message)
-                   CALL update_output(output_var, output_arr, iter, output_file)
+                   IF (.NOT. output_var % is_internal .OR. debugging_output)  THEN
+                      message = "Writing variable '" // TRIM(output_var%name) // "'."
+                      CALL report('main loop', message)
+                      CALL update_output(output_var, output_arr, iter, output_file)
  CALL run_control('time', 'write')
-
-                   DEALLOCATE(output_arr)
+                   END IF
+
+                   IF (ALLOCATED(output_arr))  DEALLOCATE(output_arr)
  CALL run_control('time', 'alloc')
 
@@ -267,4 +431,18 @@
 
              END DO ! ouput variables
+
+             IF ( group % kind == 'thermodynamics' )  THEN
+                DEALLOCATE( rho_centre )
+                IF ( .NOT. cfg % ug_is_set )  THEN
+                   DEALLOCATE( rho_north )
+                   DEALLOCATE( rho_south )
+                   DEALLOCATE( rho_east )
+                   DEALLOCATE( rho_west )
+                   DEALLOCATE( p_north )
+                   DEALLOCATE( p_south )
+                   DEALLOCATE( p_east )
+                   DEALLOCATE( p_west )
+                END IF
+             END IF
 
              IF ( group % kind == 'running average' .OR. &
@@ -273,5 +451,5 @@
                 ! accumulated COSMO-DE quantities (precipitation).
              ELSE
-                CALL report('main loop', 'Deallocating input buffer')
+                CALL report('main loop', 'Deallocating input buffer', cfg % debug)
                 DEALLOCATE(input_buffer)
              END IF
@@ -281,5 +459,5 @@
 
           IF (ALLOCATED(input_buffer))  THEN
-             CALL report('main loop', 'Deallocating input buffer')
+             CALL report('main loop', 'Deallocating input buffer', cfg % debug)
              DEALLOCATE(input_buffer)
           END IF
@@ -294,5 +472,5 @@
           END IF
 
-          CALL report('main loop', message)
+          CALL report('main loop', message, cfg % debug)
 
        END IF ! IO group % to_be_processed
Index: /palm/trunk/UTIL/inifor/src/io.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/io.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/io.f90	(revision 3395)
@@ -26,5 +26,13 @@
 ! -----------------
 ! $Id$
-! Removed unnecessary check for output file
+! Added command-line options for configuring the computation of geostrophic
+!     winds (--averagin-mode, --averaging-angle)
+! Added command-line option --input-prefix for setting input file prefixes all
+!     at once
+! Added --debug option for more verbose terminal output
+! Added option-specific *_is_set LOGICALs to indicate invocation from the
+!     command-line 
+! Improved error messages in case of empty file-name strings
+! Improved routine naming
 ! 
 ! 3183 2018-07-27 14:25:55Z suehring
@@ -148,4 +156,13 @@
        INTEGER                            ::  arg_count, i
 
+       cfg % p0_is_set = .FALSE.
+       cfg % ug_is_set = .FALSE.
+       cfg % vg_is_set = .FALSE.
+       cfg % flow_prefix_is_set = .FALSE.
+       cfg % input_prefix_is_set = .FALSE.
+       cfg % radiation_prefix_is_set = .FALSE.
+       cfg % soil_prefix_is_set = .FALSE.
+       cfg % soilmoisture_prefix_is_set = .FALSE.
+
        arg_count = COMMAND_ARGUMENT_COUNT()
        IF (arg_count .GT. 0)  THEN
@@ -170,7 +187,14 @@
              SELECT CASE( TRIM(option) )
 
+             CASE( '--averaging-mode' )
+                CALL get_option_argument( i, arg )
+                cfg % averaging_mode = TRIM(arg)
+
              CASE( '-date', '-d', '--date' )
                 CALL get_option_argument( i, arg )
                 cfg % start_date = TRIM(arg)
+
+             CASE( '--debug' )
+                cfg % debug = .TRUE.
 
              ! Elevation of the PALM-4U domain above sea level
@@ -181,4 +205,5 @@
              ! surface pressure, at z0
              CASE( '-p0', '-r', '--surface-pressure' )
+                cfg % p0_is_set = .TRUE.
                 CALL get_option_argument( i, arg )
                 READ(arg, *) cfg % p0
@@ -186,4 +211,5 @@
              ! geostrophic wind in x direction
              CASE( '-ug', '-u', '--geostrophic-u' )
+                cfg % ug_is_set = .TRUE.
                 CALL get_option_argument( i, arg )
                 READ(arg, *) cfg % ug
@@ -191,4 +217,5 @@
              ! geostrophic wind in y direction
              CASE( '-vg', '-v', '--geostrophic-v' )
+                cfg % vg_is_set = .TRUE.
                 CALL get_option_argument( i, arg )
                 READ(arg, *) cfg % vg
@@ -208,29 +235,42 @@
              CASE( '-hhl', '-l', '--hhl-file' )
                 CALL get_option_argument( i, arg )
-                 cfg % hhl_file = TRIM(arg)
+                cfg % hhl_file = TRIM(arg)
+
+             CASE( '--input-prefix')
+                CALL get_option_argument( i, arg )
+                cfg % input_prefix = TRIM(arg)
+                cfg % input_prefix_is_set = .TRUE.
+   
+             CASE( '-a', '--averaging-angle' )
+                CALL get_option_argument( i, arg )
+                READ(arg, *) cfg % averaging_angle
 
              CASE( '-static', '-t', '--static-driver' )
                 CALL get_option_argument( i, arg )
-                 cfg % static_driver_file = TRIM(arg)
+                cfg % static_driver_file = TRIM(arg)
 
              CASE( '-soil', '-s', '--soil-file')
                 CALL get_option_argument( i, arg )
-                 cfg % soiltyp_file = TRIM(arg)
+                cfg % soiltyp_file = TRIM(arg)
 
              CASE( '--flow-prefix')
                 CALL get_option_argument( i, arg )
-                 cfg % flow_prefix = TRIM(arg)
-
+                cfg % flow_prefix = TRIM(arg)
+                cfg % flow_prefix_is_set = .TRUE.
+   
              CASE( '--radiation-prefix')
                 CALL get_option_argument( i, arg )
-                 cfg % radiation_prefix = TRIM(arg)
-
+                cfg % radiation_prefix = TRIM(arg)
+                cfg % radiation_prefix_is_set = .TRUE.
+   
              CASE( '--soil-prefix')
                 CALL get_option_argument( i, arg )
-                 cfg % soil_prefix = TRIM(arg)
-
+                cfg % soil_prefix = TRIM(arg)
+                cfg % soil_prefix_is_set = .TRUE.
+   
              CASE( '--soilmoisture-prefix')
                 CALL get_option_argument( i, arg )
-                 cfg % soilmoisture_prefix = TRIM(arg)
+                cfg % soilmoisture_prefix = TRIM(arg)
+                cfg % soilmoisture_prefix_is_set = .TRUE.
 
              CASE( '-o', '--output' )
@@ -297,11 +337,11 @@
 
       all_files_present = .TRUE.
-      all_files_present = all_files_present .AND. file_present(cfg % hhl_file)
-      all_files_present = all_files_present .AND. file_present(cfg % namelist_file)
-      all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file)
+      all_files_present = all_files_present .AND. file_present(cfg % hhl_file, 'HHL')
+      all_files_present = all_files_present .AND. file_present(cfg % namelist_file, 'NAMELIST')
+      all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file, 'SOILTYP')
 
       ! Only check optional static driver file name, if it has been given.
       IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
-         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file)
+         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver')
       END IF
 
@@ -335,16 +375,44 @@
       END SELECT
 
+      SELECT CASE( TRIM(cfg % averaging_mode) )
+      CASE( 'level', 'height')
+      CASE DEFAULT
+         message = "Averaging mode '" // TRIM(cfg % averaging_mode) //&
+                   "' is not supported. " //&
+                   "Please select either 'height' or 'level', " //&
+                   "or omit the --averaging-mode option entirely, which corresponds "//&
+                   "to the latter."
+         CALL abort( 'validate_config', message ) 
+      END SELECT
+
+      IF ( cfg % ug_is_set .NEQV. cfg % vg_is_set )  THEN
+         message = "You specified only one component of the geostrophic " // &
+                   "wind. Please specify either both or none."
+         CALL abort( 'validate_config', message )
+      END IF
 
    END SUBROUTINE validate_config
 
 
-   LOGICAL FUNCTION file_present(filename)
+   LOGICAL FUNCTION file_present(filename, kind)
       CHARACTER(LEN=PATH), INTENT(IN) ::  filename
-
-      INQUIRE(FILE=filename, EXIST=file_present)
-
-      IF (.NOT. file_present)  THEN
-         message = "The given file '" // "' does not exist."
+      CHARACTER(LEN=*), INTENT(IN)    ::  kind
+
+      IF (LEN(TRIM(filename))==0)  THEN
+
+         file_present = .FALSE.
+         message = "No name was given for the " // TRIM(kind) // " file."
          CALL report('file_present', message)
+
+      ELSE
+
+         INQUIRE(FILE=filename, EXIST=file_present)
+
+         IF (.NOT. file_present)  THEN
+            message = "The given " // TRIM(kind) // " file '" //               &
+                      TRIM(filename) // "' does not exist."
+            CALL report('file_present', message)
+         END IF
+
       END IF
 
@@ -388,7 +456,17 @@
 #endif
 
-!
-!------------------------------------------------------------------------------
-!- Section 1: Write global NetCDF attributes
+!------------------------------------------------------------------------------
+!- Section 1: Define NetCDF dimensions and coordinates
+!------------------------------------------------------------------------------
+       nt = SIZE(output_file % time)
+       nx = palm_grid % nx
+       ny = palm_grid % ny
+       nz = palm_grid % nz
+       z0 = palm_grid % z0
+
+
+!
+!------------------------------------------------------------------------------
+!- Section 2: Write global NetCDF attributes
 !------------------------------------------------------------------------------
        CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
@@ -406,17 +484,12 @@
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
+       ! FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
+       ! FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
+       ! FIXME: Standard v1.9., origin_z)
+       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
 
 !
-!------------------------------------------------------------------------------
-!- Section 2: Define NetCDF dimensions and coordinates
-!------------------------------------------------------------------------------
-       nt = SIZE(output_file % time)
-       nx = palm_grid % nx
-       ny = palm_grid % ny
-       nz = palm_grid % nz
-       z0 = palm_grid % z0
-
 !
 !------------------------------------------------------------------------------
@@ -618,9 +691,16 @@
        ELSE
 
-          nbuffers = SIZE( group % in_var_list )
+          ! Allocate one input buffer per input_variable. If more quantities
+          ! have to be computed than input variables exist in this group,
+          ! allocate more buffers. For instance, in the thermodynamics group,
+          ! there are three input variabels (PP, T, Qv) and four quantities
+          ! necessart (P, Theta, Rho, qv) for the corresponding output fields
+          ! (p0, Theta, qv, ug, and vg)
+          nbuffers = MAX( group % n_inputs, group % n_output_quantities )
           ALLOCATE( buffer(nbuffers) )
  CALL run_control('time', 'alloc')
           
-          DO ivar = 1, nbuffers
+          ! Read in all input variables, leave extra buffers-if any-untouched.
+          DO ivar = 1, group % n_inputs
 
              input_var => group % in_var_list(ivar)
@@ -628,5 +708,5 @@
              ! Check wheather P or PP is present in input file
              IF (input_var % name == 'P')  THEN
-                input_var % name = TRIM( get_pressure_var(input_file) )
+                input_var % name = TRIM( get_pressure_varname(input_file) )
  CALL run_control('time', 'read')
              END IF
@@ -668,5 +748,5 @@
 !> perturbation, 'PP', and returns the appropriate string.
 !------------------------------------------------------------------------------!
-    CHARACTER(LEN=2) FUNCTION get_pressure_var(input_file) RESULT(var)
+    CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
        CHARACTER(LEN=*) ::  input_file
        INTEGER          ::  ncid, varid
@@ -692,5 +772,5 @@
        CALL check(nf90_close(ncid))
 
-    END FUNCTION get_pressure_var
+    END FUNCTION get_pressure_varname
 
 
@@ -805,5 +885,5 @@
                                    start=start(1:ndim+1) ) )
 
-       CASE ( 'constant scalar profile' )
+       CASE ( 'constant scalar profile', 'geostrophic', 'internal profile' )
 
           CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
Index: /palm/trunk/UTIL/inifor/src/transform.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/transform.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/transform.f90	(revision 3395)
@@ -26,4 +26,11 @@
 ! -----------------
 ! $Id$
+! Switched addressing of averaging regions from index bounds to list of columns
+! Added routines for the computation of geostrophic winds including:
+!    - the downward extrapolation of density (linear) and
+!    - pressure (hydrostatic equation) and
+!    - rotation of geostrophic wind components to PALM frame of reference
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Introduced new PALM grid stretching
 ! Removed unnecessary subroutine parameters
@@ -51,5 +58,5 @@
     USE control
     USE defs,                                                                  &
-        ONLY: TO_DEGREES, TO_RADIANS, PI, dp
+        ONLY: G, TO_DEGREES, TO_RADIANS, PI, dp
     USE types
     USE util,                                                                  &       
@@ -176,30 +183,32 @@
 
 
-    SUBROUTINE average_2d(in_arr, out_arr, imin, imax, jmin, jmax)
-       REAL(dp), INTENT(IN)  ::  in_arr(0:,0:,0:)
-       REAL(dp), INTENT(OUT) ::  out_arr(0:)
-       INTEGER, INTENT(IN)   ::  imin, imax, jmin, jmax
-
-       INTEGER  ::  i, j, k
+    SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
+       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
+       REAL(dp), INTENT(OUT)             ::  out_arr(0:)
+       INTEGER, INTENT(IN), DIMENSION(:) ::  ii, jj
+
+       INTEGER  ::  i, j, k, l
        REAL(dp) ::  ni
-       
-       IF (imin < 0)  CALL abort('average_2d', "imin < 0.")
-       IF (jmin < 0)  CALL abort('average_2d', "jmin < 0.")
-       IF (imax > UBOUND(in_arr, 1))  CALL abort('average_2d', "imax out of i bound.")
-       IF (jmax > UBOUND(in_arr, 2))  CALL abort('average_2d', "jmax out of j bound.")
+
+       IF (SIZE(ii) .NE. SIZE(jj))  THEN
+          message = "Length of 'ii' and 'jj' index lists do not match." //     &
+             NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " //        &
+             NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
+          CALL abort('average_2d', message)
+       END IF
 
        DO k = 0, UBOUND(out_arr, 1)
 
           out_arr(k) = 0.0_dp
-          DO j = jmin, jmax
-          DO i = imin, imax
-             out_arr(k) = out_arr(k) + in_arr(i, j, k)
+          DO l = 1, UBOUND(ii, 1)
+             i = ii(l)
+             j = jj(l)
+             out_arr(k) = out_arr(k) +&
+                         in_arr(i, j, k)
           END DO
-          END DO
-
-       END DO
-   
-       ! devide by number of grid points
-       ni = 1.0_dp / ( (imax - imin + 1) * (jmax - jmin + 1) )
+
+       END DO
+
+       ni = 1.0_dp / SIZE(ii)
        out_arr(:) = out_arr(:) * ni
 
@@ -212,14 +221,14 @@
        REAL(dp), DIMENSION(:,:,:), INTENT(OUT) ::  palm_array
        REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::  intermediate_array
-       INTEGER ::  nx, ny, nz
+       INTEGER ::  nx, ny, nlev
 
        nx = palm_intermediate % nx
        ny = palm_intermediate % ny
-       nz = palm_intermediate % nz ! nlev
+       nlev = palm_intermediate % nz ! nlev
 
        ! Interpolate from COSMO-DE to intermediate grid. Allocating with one
        ! less point in the vertical, since scalars like T have 50 instead of 51
        ! points in COSMO-DE.
-       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1)) !
+       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) !
 
        CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
@@ -234,29 +243,144 @@
 
 
-    SUBROUTINE average_profile(source_array, profile_array, imin, imax, jmin, jmax,&
-                               palm_intermediate, palm_grid)
-       TYPE(grid_definition), INTENT(IN)          ::  palm_intermediate, palm_grid
+    SUBROUTINE average_profile(source_array, profile_array, avg_grid)
+       TYPE(grid_definition), INTENT(IN)          ::  avg_grid
        REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array
-       INTEGER, INTENT(IN)                        ::  imin, imax, jmin, jmax
-       REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) ::  profile_array
-       REAL(dp), DIMENSION(:,:,:), ALLOCATABLE    ::  intermediate_array
-       INTEGER                                    ::  nx, ny, nz
-
-       nx = palm_intermediate % nx
-       ny = palm_intermediate % ny
-       nz = palm_intermediate % nz
-       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1))
-       intermediate_array(:,:,:) = 0.0_dp
-
-       ! average input array to intermediate profile
-       CALL average_2d(source_array, intermediate_array(0,0,:), imin, imax, jmin, jmax)
-
-       ! vertically interpolate to ouput array
-       CALL interpolate_1d(intermediate_array, profile_array, palm_grid)
-
-       DEALLOCATE(intermediate_array)
+       REAL(dp), DIMENSION(:), INTENT(OUT)        ::  profile_array
+
+       INTEGER ::  i_source, j_source, k_profile, k_source, l, m
+
+       REAL ::  ni_columns
+
+       profile_array(:) = 0.0_dp
+
+       DO l = 1, avg_grid % n_columns
+          i_source = avg_grid % iii(l)
+          j_source = avg_grid % jjj(l)
+
+          DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) ! PALM levels
+
+             DO m = 1, 2 ! vertical interpolation neighbours
+
+                k_source = avg_grid % kkk(l, k_profile, m)
+
+                profile_array(k_profile) = profile_array(k_profile)    &
+                   + avg_grid % w(l, k_profile, m)                             &
+                   * source_array(i_source, j_source, k_source)
+
+             END DO ! m, vertical interpolation neighbours
+
+          END DO    ! k_profile, PALM levels
+
+       END DO       ! l, horizontal neighbours
+
+       ni_columns = 1.0_dp / avg_grid % n_columns
+       profile_array(:) = profile_array(:) * ni_columns
+
+       ! Extrapolate constant to the bottom
+       profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
 
     END SUBROUTINE average_profile
 
+
+    SUBROUTINE extrapolate_density(rho, avg_grid)
+       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  rho
+       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
+
+       REAL(dp) ::  drhodz, dz, zk, rhok
+       INTEGER  ::  k_min
+
+       k_min  = avg_grid % k_min
+       zk     = avg_grid % z(k_min)
+       rhok   = rho(k_min)
+       dz     = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
+       drhodz = (rho(k_min + 1) - rho(k_min)) / dz
+
+       rho(1:k_min-1) = rhok + drhodz * (avg_grid % z(1:k_min-1) - zk)
+
+    END SUBROUTINE extrapolate_density
+
+
+    SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
+       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
+       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
+       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
+
+       REAL(dp) ::  drhodz, dz, zk, rhok
+       INTEGER  ::  k, k_min
+
+       k_min = avg_grid % k_min
+       zk    = avg_grid % z(k_min)
+       rhok  = rho(k_min)
+       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
+       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
+
+       DO k = 1, k_min-1
+          p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz,         &
+                                           avg_grid % z(k), G)
+       END DO
+
+    END SUBROUTINE extrapolate_pressure
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Takes the averaged pressure profile <p> and sets the lowest entry to the
+!> extrapolated pressure at the surface.
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_surface_pressure(p, rho, avg_grid)
+       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
+       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
+       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
+
+       REAL(dp) ::  drhodz, dz, zk, rhok
+       INTEGER  ::  k, k_min
+
+       k_min = avg_grid % k_min
+       zk    = avg_grid % z(k_min)
+       rhok  = rho(k_min)
+       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
+       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
+
+       p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz,            &
+                                        0.0, G)
+
+    END SUBROUTINE get_surface_pressure
+
+
+    FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
+
+       REAL(dp), INTENT(IN)  ::  pk, zk, rhok, drhodz, g, z
+       REAL(dp) ::  p
+
+       p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) )
+
+    END FUNCTION constant_density_pressure
+
+!-----------------------------------------------------------------------------!
+! Description:
+! -----------
+!> This routine computes profiles of the two geostrophic wind components ug and
+!> vg.
+!-----------------------------------------------------------------------------!
+    SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3,    &
+                                 Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg)
+
+    REAL(dp), DIMENSION(:), INTENT(IN)  ::  p_north, p_south, p_east, p_west,  &
+                                            rho
+    REAL(dp), INTENT(IN)                ::  f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g
+    REAL(dp), DIMENSION(:), INTENT(OUT) ::  ug, vg
+    REAL(dp)                            ::  facx, facy
+
+    facx = 1.0_dp / (Lx * f3)
+    facy = 1.0_dp / (Ly * f3)
+    ug(:) = - facy / rho(:) * (p_north(:) - p_south(:))
+    vg(:) =   facx / rho(:) * (p_east(:) - p_west(:))
+
+    CALL rotate_vector_field(                                                  &
+       ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)&
+    )
+
+    END SUBROUTINE geostrophic_winds
 
 
@@ -419,4 +543,54 @@
        
 
+    SUBROUTINE rotate_vector_field(x, y, angle)
+       REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y  !< x and y coodrinate in arbitrary units
+       REAL(dp), INTENT(IN)                  :: angle !< rotation angle [deg]
+
+       INTEGER  :: i
+       REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2)
+
+       sine = SIN(angle * TO_RADIANS)
+       cosine = COS(angle * TO_RADIANS)
+       ! RESAHPE() fills columns first, so the rotation matrix becomes
+       !
+       ! rotation = [ cosine   -sine  ]
+       !            [  sine    cosine ]
+       rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
+
+       DO i = LBOUND(x, 1), UBOUND(x, 1)
+
+          v_rot(:) = MATMUL(rotation, (/x(i), y(i)/))
+
+          x(i) = v_rot(1)
+          y(i) = v_rot(2)
+
+       END DO
+
+    END SUBROUTINE rotate_vector_field
+
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine computes the local meridian convergence between a rotated-pole
+!> and a geographical system using the Eq. (6) given in the DWD manual
+!> 
+!>    Baldauf et al. (2018), Beschreibung des operationelle KuÌrzestfrist-
+!>    vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die
+!>    Datenbanken des DWD.
+!>    https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2
+!>
+    FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)          &
+       RESULT(delta)
+
+       REAL(dp), INTENT(IN) ::  phi_n, lam_n, phi_g, lam_g
+       REAL(dp)             ::  delta
+
+       delta = atan2( COS(phi_n) * SIN(lam_n - lam_g),                         &
+                      COS(phi_g) * SIN(phi_n) -                                &
+                      SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) )
+
+    END FUNCTION meridian_convergence_rotated
 
 !------------------------------------------------------------------------------!
@@ -512,5 +686,6 @@
 
     
-    SUBROUTINE find_vertical_neighbours_and_weights(palm_grid, palm_intermediate)
+    SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
+                                                            palm_intermediate )
        TYPE(grid_definition), INTENT(INOUT) ::  palm_grid
        TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
@@ -528,6 +703,6 @@
 
        ! in each column of the fine grid, find vertical neighbours of every cell
+       DO j = 0, ny
        DO i = 0, nx
-       DO j = 0, ny
 
           k_intermediate = 0
@@ -536,5 +711,5 @@
           column_top  = palm_intermediate % h(i,j,nlev)
 
-          ! scan through palm_grid column until and set neighbour indices in
+          ! scan through palm_grid column and set neighbour indices in
           ! case current_height is either below column_base, in the current
           ! cell, or above column_top. Keep increasing current cell index until
@@ -580,4 +755,101 @@
                    h_top    = palm_intermediate % h(i,j,k_intermediate+1)
                    h_bottom = palm_intermediate % h(i,j,k_intermediate)
+                   point_is_in_current_cell = (                                &
+                      current_height >= h_bottom .AND.                         &
+                      current_height <  h_top                                  &
+                   )
+                END DO
+
+                IF (k_intermediate > nlev-1)  THEN
+                   message = "Index " // TRIM(str(k_intermediate)) //          &
+                             " is above intermediate grid range."
+                   CALL abort('find_vertical_neighbours', message)
+                END IF
+   
+                palm_grid % kk(i,j,k,1) = k_intermediate
+                palm_grid % kk(i,j,k,2) = k_intermediate + 1
+
+                ! copmute vertical weights
+                weight = (h_top - current_height) / (h_top - h_bottom)
+                palm_grid % w_verti(i,j,k,1) = weight
+                palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
+             END IF
+
+          END DO
+
+       END DO
+       END DO
+
+    END SUBROUTINE find_vertical_neighbours_and_weights_interp
+
+
+    SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid )
+       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
+
+       INTEGER  ::  i, j, k_palm, k_intermediate, l, nlev
+       LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
+                    point_is_in_current_cell
+       REAL(dp) ::  current_height, column_base, column_top, h_top, h_bottom,  &
+                    weight
+
+
+       avg_grid % k_min = LBOUND(avg_grid % z, 1)
+
+       nlev = SIZE(avg_grid % cosmo_h, 3)
+
+       ! in each column of the fine grid, find vertical neighbours of every cell
+       DO l = 1, avg_grid % n_columns
+
+          i = avg_grid % iii(l)
+          j = avg_grid % jjj(l)
+
+          column_base = avg_grid % cosmo_h(i,j,1)
+          column_top  = avg_grid % cosmo_h(i,j,nlev)
+
+          ! scan through avg_grid column until and set neighbour indices in
+          ! case current_height is either below column_base, in the current
+          ! cell, or above column_top. Keep increasing current cell index until
+          ! the current cell overlaps with the current_height.
+          k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based.
+          DO k_palm = 1, avg_grid % nz
+
+             ! Memorize the top and bottom boundaries of the coarse cell and the
+             ! current height within it
+             current_height = avg_grid % z(k_palm) + avg_grid % z0
+             h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
+             h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
+
+             point_is_above_grid = (current_height > column_top) !22000m, very unlikely
+             point_is_below_grid = (current_height < column_base)
+
+             point_is_in_current_cell = (                                      &
+                current_height >= h_bottom .AND.                               &
+                current_height <  h_top                                        &
+             )
+
+             ! set default weights
+             avg_grid % w(l,k_palm,1:2) = 0.0_dp
+
+             IF (point_is_above_grid)  THEN
+
+                avg_grid % kkk(l,k_palm,1:2) = nlev
+                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
+
+                message = "PALM-4U grid extends above COSMO-DE model top."
+                CALL abort('find_vertical_neighbours_and_weights_average', message)
+
+             ELSE IF (point_is_below_grid)  THEN
+
+                avg_grid % kkk(l,k_palm,1:2) = 0
+                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
+                avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min)
+             ELSE
+                ! cycle through intermediate levels until current
+                ! intermediate-grid cell overlaps with current_height
+                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
+                   k_intermediate = k_intermediate + 1
+
+                   h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
+                   h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
                    point_is_in_current_cell = (                                &
                       current_height >= h_bottom .AND.                         &
@@ -594,19 +866,17 @@
                 END IF
    
-                palm_grid % kk(i,j,k,1) = k_intermediate
-                palm_grid % kk(i,j,k,2) = k_intermediate + 1
+                avg_grid % kkk(l,k_palm,1) = k_intermediate
+                avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
 
                 ! copmute vertical weights
                 weight = (h_top - current_height) / (h_top - h_bottom)
-                palm_grid % w_verti(i,j,k,1) = weight
-                palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
+                avg_grid % w(l,k_palm,1) = weight
+                avg_grid % w(l,k_palm,2) = 1.0_dp - weight
              END IF
 
-          END DO
-
-       END DO
-       END DO
-
-    END SUBROUTINE find_vertical_neighbours_and_weights
+          END DO ! k, PALM levels
+       END DO ! l, averaging columns
+ 
+    END SUBROUTINE find_vertical_neighbours_and_weights_average
 
 !------------------------------------------------------------------------------!
Index: /palm/trunk/UTIL/inifor/src/types.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/types.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/types.f90	(revision 3395)
@@ -26,4 +26,11 @@
 ! -----------------
 ! $Id$
+! Added *_is_set LOGICALs to inifor_config type to indicate option invocation
+!     from the command-line 
+! Added 1D index vertical weights lists to support addressing averaging regions
+!     by list of columns instead of index bounds
+! 
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Introduced new PALM grid stretching:
 ! - Converted vertical grid_definition coordinte variables to pointers
@@ -66,8 +73,10 @@
 
     CHARACTER(LEN=SNAME) ::  flow_prefix          !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  input_prefix         !< Prefix of all input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses
     CHARACTER(LEN=SNAME) ::  soil_prefix          !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
-    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses
     CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses
 
+    CHARACTER(LEN=SNAME) ::  averaging_mode
     CHARACTER(LEN=SNAME) ::  bc_mode
     CHARACTER(LEN=SNAME) ::  ic_mode
@@ -77,15 +86,34 @@
     REAL(dp)             ::  ug
     REAL(dp)             ::  vg
-    REAL(dp)             ::  z0 !< Elevation of the PALM-4U domain above sea level [m]
+    REAL(dp)             ::  z0                   !< Elevation of the PALM-4U domain above sea level [m]
+    REAL(dp)             ::  averaging_angle      !< latitudal and longitudal width of averaging regions [deg]
+    
+
+    LOGICAL              ::  debug
+    LOGICAL              ::  p0_is_set
+    LOGICAL              ::  ug_is_set
+    LOGICAL              ::  vg_is_set
+    LOGICAL              ::  flow_prefix_is_set          !< 
+    LOGICAL              ::  input_prefix_is_set         !< 
+    LOGICAL              ::  radiation_prefix_is_set     !< 
+    LOGICAL              ::  soil_prefix_is_set          !< 
+    LOGICAL              ::  soilmoisture_prefix_is_set  !< 
  END TYPE inifor_config
 
  TYPE grid_definition
     CHARACTER(LEN=SNAME)  ::  name(3)       !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/)
+    CHARACTER(LEN=SNAME)  ::  kind          !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/)
+    INTEGER               ::  k_min         !< Index of lowest PALM grid level that is not cut by local COSMO orography; vertically separates interpolation and extrapolation region.
     INTEGER               ::  nx            !< number of gridpoints in the first dimension
     INTEGER               ::  ny            !< number of gridpoints in the second dimension
-    INTEGER               ::  nz            !< number of gridpoints in the third dimension
+    INTEGER               ::  nz            !< number of gridpoints in the third dimension, used for PALM points
+    INTEGER               ::  nlev          !< number of COSMO grid levels
+    INTEGER               ::  n_columns     !< number of averaging columns of the source grid 
     INTEGER, ALLOCATABLE  ::  ii(:,:,:)     !< Given a point (i,j,k) in the PALM-4U grid, ii(i,j,l) gives the x index of the l'th horizontl neighbour on the COSMO-DE grid.
     INTEGER, ALLOCATABLE  ::  jj(:,:,:)     !< Given a point (i,j,k) in the PALM-4U grid, jj(i,j,l) gives the y index of the l'th horizontl neighbour on the COSMO-DE grid.
     INTEGER, ALLOCATABLE  ::  kk(:,:,:,:)   !< Given a point (i,j,k) in the PALM-4U grid, kk(i,j,k,l) gives the z index of the l'th vertical neighbour in the intermediate grid.
+    INTEGER, ALLOCATABLE  ::  iii(:)        !< profile averaging neighbour indices 
+    INTEGER, ALLOCATABLE  ::  jjj(:)        !< profile averaging neighbour indices 
+    INTEGER, ALLOCATABLE  ::  kkk(:,:,:)    !< indices of vertical interpolation neightbours, kkk(<source column>, <PALM k level>, <neighbour index>)
     REAL(dp)              ::  lx            !< domain length in the first dimension [m]
     REAL(dp)              ::  ly            !< domain length in the second dimension [m]
@@ -97,4 +125,5 @@
     REAL(dp), POINTER     ::  z(:)          !< coordinates of cell centers in z direction [m]
     REAL(dp), ALLOCATABLE ::  h(:,:,:)      !< heights grid point for intermediate grids [m]
+    REAL(dp), POINTER     ::  cosmo_h(:,:,:)!< pointer to appropriate COSMO level heights (scalar/w) [m]
     REAL(dp), POINTER     ::  hhl(:,:,:)    !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from 
     REAL(dp), POINTER     ::  hfl(:,:,:)    !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
@@ -115,4 +144,5 @@
     REAL(dp), ALLOCATABLE ::  w_horiz(:,:,:)   !< weights for bilinear horizontal interpolation
     REAL(dp), ALLOCATABLE ::  w_verti(:,:,:,:) !< weights for linear vertical interpolation
+    REAL(dp), ALLOCATABLE ::  w(:,:,:)      !< vertical interpolation weights, w(<source_column>, <PALM k level>, <neighbour index>) [-]
  END TYPE grid_definition
 
@@ -150,8 +180,10 @@
     CHARACTER(LEN=SNAME)                  ::  task                      !< Processing task that generates this variable, e.g. 'interpolate_2d' or 'average profile'
     LOGICAL                               ::  to_be_processed = .FALSE. !< INIFOR flag indicating whether variable shall be processed
+    LOGICAL                               ::  is_internal = .FALSE.     !< INIFOR flag indicating whether variable shall be written to netCDF file (.FALSE.) or kept for later (.TRUE.)
     LOGICAL                               ::  is_read = .FALSE.         !< INIFOR flag indicating whether variable has been read
-    LOGICAL                               ::  is_upside_down  = .FALSE. !< INIFOR flag indicating whether variable shall be processed
+    LOGICAL                               ::  is_upside_down  = .FALSE. !< INIFOR flag indicating whether vertical dimension is reversed (typically the case with COSMO-DE atmospheric fields)
     TYPE(grid_definition), POINTER        ::  grid                      !< Pointer to the corresponding output grid
     TYPE(grid_definition), POINTER        ::  intermediate_grid         !< Pointer to the corresponding intermediate grid
+    TYPE(grid_definition), POINTER        ::  averaging_grid         !< Pointer to the corresponding intermediate grid
  END TYPE nc_var
 
@@ -159,5 +191,7 @@
  TYPE io_group                                          !< Input/Output group, groups together output variabels that share their input variables. For instance, all boundary surfaces and initialization fields of the potential temperature are base on T and p.
     INTEGER                          ::  nt             !< maximum number of output time steps across all output variables
-    INTEGER                          ::  nv             !< number of output variables
+    INTEGER                          ::  nv             !< number of netCDF output variables
+    INTEGER                          ::  n_inputs       !< number of input variables
+    INTEGER                          ::  n_output_quantities !< number of physical quantities required for computing netCDF output variables
     CHARACTER(LEN=SNAME)             ::  kind           !< kind of I/O group
     CHARACTER(LEN=PATH), ALLOCATABLE ::  in_files(:)    !< list of nt input files
Index: /palm/trunk/UTIL/inifor/src/util.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/util.f90	(revision 3394)
+++ /palm/trunk/UTIL/inifor/src/util.f90	(revision 3395)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! New routines for computing potential temperature and moist air density
+! Increased number of digits in real-to-str conversion
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Improved real-to-string conversion
 ! 
@@ -48,4 +52,6 @@
     USE defs,                                                                  &
         ONLY :  dp, PI, DATE, SNAME
+    USE types,                                                                 &
+        ONLY :  grid_definition
 
     IMPLICIT NONE
@@ -281,4 +287,41 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Converts the absolute temperature to the potential temperature in place using
+!> the identity a^b = e^(b ln(a)).
+!>
+!>     theta = T * (p_ref/p)^(R/c_p) = T * e^( R/c_p * ln(p_ref/p) )
+!------------------------------------------------------------------------------!
+    SUBROUTINE potential_temperature(t, p, p_ref, r, cp)
+       REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) ::  t
+       REAL(dp), DIMENSION(:,:,:), INTENT(IN)    ::  p
+       REAL(dp), INTENT(IN)                      ::  p_ref, r, cp
+       REAL(dp)                                  ::  rcp
+
+       rcp = r/cp
+       t(:,:,:) =  t(:,:,:) * EXP( rcp * LOG(p_ref / p(:,:,:)) )
+
+    END SUBROUTINE potential_temperature
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the density in place of the given temperature (t_rho).
+!------------------------------------------------------------------------------!
+   SUBROUTINE moist_density(t_rho, p, qv, rd, rv)
+       REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) ::  t_rho
+       REAL(dp), DIMENSION(:,:,:), INTENT(IN)    ::  p, qv
+       REAL(dp), INTENT(IN)                      ::  rd, rv
+
+       t_rho(:,:,:) = p(:,:,:) / (                                             &
+          (rv * qv(:,:,:) + rd * (1.0_dp - qv(:,:,:))) * t_rho(:,:,:)          &
+       )
+
+    END SUBROUTINE moist_density
+
+
     ! Convert a real number to a string in scientific notation
     ! showing four significant digits.
@@ -298,9 +341,9 @@
 
 
-    CHARACTER(LEN=12) FUNCTION real_to_str_f(val)
+    CHARACTER(LEN=16) FUNCTION real_to_str_f(val)
 
         REAL(dp), INTENT(IN) ::  val
 
-        WRITE(real_to_str_f, '(F12.4)') val
+        WRITE(real_to_str_f, '(F16.8)') val
         real_to_str_f = ADJUSTL(real_to_str_f)
 
Index: lm/trunk/UTIL/inifor/tests/test-stretching.f90
===================================================================
--- /palm/trunk/UTIL/inifor/tests/test-stretching.f90	(revision 3394)
+++ 	(revision )
@@ -1,100 +1,0 @@
-!> @file tests/test-stretching.f90
-!------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017-2018 Leibniz Universitaet Hannover
-! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
-!------------------------------------------------------------------------------!
-!
-! Current revisions:
-! -----------------
-! 
-! 
-! Former revisions:
-! -----------------
-! $Id$
-! Initial revision
-! 
-! 
-! 
-! 
-!
-! Authors:
-! --------
-! @author Eckhard Kadasch
-!
-! Description:
-! ------------
-!> This program tests INIFOR's implementation of PALM's grid stretching.
-!> 
-!------------------------------------------------------------------------------!
- PROGRAM test_stretching
-
-    USE defs, &
-        ONLY :  dp
-
-    USE grid, &
-        ONLY :  stretched_z
-
-    USE test_utils
-    
-    IMPLICIT NONE
-
-    CHARACTER(LEN=*), PARAMETER ::  title = "stretched grid"
-    LOGICAL                     ::  res
-
-    INTEGER, PARAMETER ::  nz = 9
-    INTEGER            ::  k  = nz - 2
-
-    REAL(dp) ::  z(1:nz)
-    REAL(dp) ::  dz(10)            = -1.0_dp
-    REAL(dp) ::  dz_max            = 1000.0_dp
-    REAL(dp) ::  dz_stretch_factor = 1.08_dp
-    REAL(dp) ::  dz_stretch_level  = 2.0_dp
-    REAL(dp) ::  dz_stretch_level_start(9) = -9999999.9_dp
-    REAL(dp) ::  dz_stretch_level_end(9) = 9999999.9_dp
-    REAL(dp) ::  dz_stretch_factor_array(9) = 1.08_dp
-
-    CALL begin_test(title, res)
-
-    ! Arange
-    z(:)   = 0.0_dp
-    dz(1)  = 1.0_dp
-
-    ! Act
-    CALL stretched_z(z, 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)
-
-    ! Assert, that the total distance covered by the stretched region 
-    ! matches the therotetial distance, i.e. the sum over the finite
-    ! exponential series
-    !
-    !          Sum_{i=0}^n( dz[i] ) = dz[0] * Sum_{i=0}^n( f^i )
-    !                               = dz[0] * (1 - f^(i+1)) / (1-f)
-    !
-    ! with f being the stretch factor.
-    res = res .AND. &
-          assert_equal( (/ z(UBOUND(z, 1)) - z(1)              /),             &
-                        (/ (1.0_dp - dz_stretch_factor**(k+1)) /               &
-                           (1.0_dp - dz_stretch_factor)        /),             &
-                        'length of stretched grid' )
-
-    CALL end_test(title, res)
-
- END PROGRAM test_stretching
