Index: palm/trunk/UTIL/inifor/src/inifor.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor.f90	(revision 3557)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+! 
+! 3537 2018-11-20 10:53:14Z eckhard
 ! Print version number on program start
 ! 
@@ -54,5 +57,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch, (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -84,26 +87,30 @@
     IMPLICIT NONE
     
-    INTEGER                                 ::  igroup
-    INTEGER                                 ::  ivar
-    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.
+    INTEGER ::  igroup !< loop index for IO groups
+    INTEGER ::  ivar   !< loop index for output variables
+    INTEGER ::  iter   !< loop index for time steps
+
+    REAL(dp), ALLOCATABLE, DIMENSION(:,:,:)     ::  output_arr !< array buffer for interpolated quantities
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre !< density profile of the centre averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_arr     !< geostrophic wind in x direction
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_arr     !< geostrophic wind in y direction
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north  !< density profile of the northern averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south  !< density profile of the southern averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east   !< density profile of the eastern averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west   !< density profile of the western averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north    !< pressure profile of the northern averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south    !< pressure profile of the southern averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east     !< pressure profile of the eastern averaging domain
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west     !< pressure profile of the western averaging domain
+
+    REAL(dp), POINTER, DIMENSION(:) ::  internal_arr !< pointer to the currently processed internal array (density, pressure)
+    REAL(dp), POINTER, DIMENSION(:) ::  ug_vg_arr    !< pointer to the currently processed geostrophic wind component
+
+    TYPE(nc_var), POINTER        ::  output_var      !< pointer to the currently processed output variable
+    TYPE(io_group), POINTER      ::  group           !< pointer to the currently processed IO group
+    TYPE(container), ALLOCATABLE ::  input_buffer(:) !< buffer of the current IO group's input arrays
+
+    LOGICAL, SAVE ::  ug_vg_have_been_computed = .FALSE. !< flag for managing geostrophic wind allocation and computation
+    LOGICAL, SAVE ::  debugging_output = .TRUE.          !< flag controllging output of internal variables
     
 !> \mainpage About INIFOR
@@ -116,22 +123,27 @@
     CALL report('main_loop', 'Running INIFOR version ' // VERSION)
 
-    ! Initialize INIFOR's parameters from command-line interface and namelists
+!
+!-- Initialize INIFOR's parameters from command-line interface and namelists
     CALL setup_parameters()
 
-    ! Initialize all grids, including interpolation neighbours and weights
+!
+!-- Initialize all grids, including interpolation neighbours and weights
     CALL setup_grids()
  CALL run_control('time', 'init')
 
-    ! Initialize the netCDF output file and define dimensions
+!
+!-- Initialize the netCDF output file and define dimensions
     CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date,    &
                                  origin_lon, origin_lat)
  CALL run_control('time', 'write')
 
-    ! Set up the tables containing the input and output variables and set
-    ! the corresponding netCDF dimensions for each output variable
+!
+!-- Set up the tables containing the input and output variables and set
+!-- the corresponding netCDF dimensions for each output variable
     CALL setup_variable_tables(cfg % ic_mode)
  CALL run_control('time', 'write')
 
-    ! Add the output variables to the netCDF output file
+!
+!-- Add the output variables to the netCDF output file
     CALL setup_netcdf_variables(output_file % name, output_var_table,          &
                                 cfg % debug)
@@ -143,5 +155,4 @@
 !- Section 2: Main loop
 !------------------------------------------------------------------------------
-
     DO igroup = 1, SIZE(io_group_list)
 
@@ -345,9 +356,9 @@
  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.
+!
+!--                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' )
 
@@ -360,9 +371,10 @@
                             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, &
+                            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
@@ -372,5 +384,6 @@
                       END IF
 
-                      ! Prepare output of geostrophic winds
+!
+!--                   Prepare output of geostrophic winds
                       SELECT CASE(TRIM(output_var % name))
                       CASE ('ls_forcing_ug')
@@ -397,10 +410,4 @@
                       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('nudging_tau')
                           output_arr(1, 1, :) = NUDGING_TAU
@@ -418,5 +425,4 @@
                                 "has not been implemented, yet."
                       CALL abort('main loop', message)
-                      !ALLOCATE( output_arr( 1, 1, 1:nz ) )
 
                    CASE DEFAULT
@@ -444,5 +450,7 @@
                 END IF
 
-             END DO ! ouput variables
+!
+!--          output variable loop
+             END DO
 
              ug_vg_have_been_computed = .FALSE.
@@ -463,8 +471,9 @@
              END IF
 
+!
+!--          Keep input buffer around for averaged (radiation) and 
+!--          accumulated COSMO quantities (precipitation).
              IF ( group % kind == 'running average' .OR. &
                   group % kind == 'accumulated' )  THEN
-                ! Keep input buffer around for averaged (radiation) and 
-                ! accumulated COSMO-DE quantities (precipitation).
              ELSE
                 CALL report('main loop', 'Deallocating input buffer', cfg % debug)
@@ -473,5 +482,7 @@
  CALL run_control('time', 'alloc')
 
-          END DO ! time steps / input files 
+!
+!--       time steps / input files loop
+          END DO
 
           IF (ALLOCATED(input_buffer))  THEN
@@ -491,7 +502,11 @@
           CALL report('main loop', message, cfg % debug)
 
-       END IF ! IO group % to_be_processed
-
-    END DO ! IO groups
+!
+!--    IO group % to_be_processed conditional
+       END IF
+
+!
+!-- IO groups loop
+    END DO
 
 !------------------------------------------------------------------------------
Index: palm/trunk/UTIL/inifor/src/inifor_control.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_control.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor_control.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+! 
+! 
+! 3447 2018-10-29 15:52:54Z eckhard
 ! Renamed source files for compatibilty with PALM build system
 !
@@ -44,5 +48,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -61,19 +65,31 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=5000) ::  message = ''
+    CHARACTER (LEN=5000) ::  message = '' !< log message buffer
 
  CONTAINS
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> 
+!> report() is INIFOR's general logging routine. It prints the given 'message'
+!> to the terminal and writes it to the INIFOR log file.
+!> 
+!> You can use this routine to log events across INIFOR's code to log. For
+!> warnings and abort messages, you may use the dedicated routines warn() and
+!> abort() in this module. Both use report() and add specific behaviour to it.
+!------------------------------------------------------------------------------!
     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
+       CHARACTER(LEN=*), INTENT(IN)  ::  routine !< name of calling subroutine of function
+       CHARACTER(LEN=*), INTENT(IN)  ::  message !< log message
+       LOGICAL, OPTIONAL, INTENT(IN) ::  debug   !< flag the current message as debugging message
+
+       INTEGER                       ::  u                     !< Fortran file unit for the log file
+       LOGICAL, SAVE                 ::  is_first_run = .TRUE. !< control flag for file opening mode
+       LOGICAL                       ::  suppress_message      !< control falg for additional debugging log
+
+
+       IF ( is_first_run )  THEN
           OPEN( NEWUNIT=u, FILE='inifor.log', STATUS='replace' )
           is_first_run = .FALSE.
@@ -84,9 +100,9 @@
 
        suppress_message = .FALSE.
-       IF (PRESENT(debug))  THEN
-          IF (.NOT. debug)  suppress_message = .TRUE.
+       IF ( PRESENT(debug) )  THEN
+          IF ( .NOT. debug )  suppress_message = .TRUE.
        END IF
 
-       IF (.NOT. suppress_message)  THEN
+       IF ( .NOT. suppress_message )  THEN
           PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
           WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
@@ -98,8 +114,18 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> 
+!> warn() prepends "WARNING:" the given 'message' and prints the result to the
+!> terminal and writes it to the INIFOR logfile.
+!> 
+!> You can use this routine for messaging issues, that still allow INIFOR to
+!> continue.
+!------------------------------------------------------------------------------!
     SUBROUTINE warn(routine, message)
 
-       CHARACTER(LEN=*), INTENT(IN) ::  routine
-       CHARACTER(LEN=*), INTENT(IN) ::  message
+       CHARACTER(LEN=*), INTENT(IN) ::  routine !< name of calling subroutine or function
+       CHARACTER(LEN=*), INTENT(IN) ::  message !< log message
 
        CALL report(routine, "WARNING: " // TRIM(message))
@@ -108,8 +134,18 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> 
+!> abort() prepends "ERROR:" the given 'message' and prints the result to the
+!> terminal, writes it to the INIFOR logfile, and exits INIFOR.
+!> 
+!> You can use this routine for messaging issues, that are critical and prevent
+!> INIFOR from continueing.
+!------------------------------------------------------------------------------!
     SUBROUTINE abort(routine, message)
 
-       CHARACTER(LEN=*), INTENT(IN) ::  routine
-       CHARACTER(LEN=*), INTENT(IN) ::  message
+       CHARACTER(LEN=*), INTENT(IN) ::  routine !< name of calling subroutine or function
+       CHARACTER(LEN=*), INTENT(IN) ::  message !< log message
 
        CALL report(routine, "ERROR: " // TRIM(message) // " Stopping.")
@@ -119,4 +155,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> 
+!> print_version() prints the INIFOR version number and copyright notice.
+!------------------------------------------------------------------------------!
     SUBROUTINE print_version()
        PRINT *, "INIFOR " // VERSION
@@ -125,15 +167,26 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> 
+!> run_control() measures the run times of various parts of INIFOR and
+!> accumulates them in timing budgets.
+!------------------------------------------------------------------------------!
     SUBROUTINE run_control(mode, budget)
 
-       CHARACTER(LEN=*), INTENT(IN) ::  mode, budget
-       REAL(dp), SAVE               ::  t0, t1
-       REAL(dp), SAVE               ::  t_comp=0.0_dp, &
-                                        t_alloc=0.0_dp, &
-                                        t_init=0.0_dp, &
-                                        t_read=0.0_dp, &
-                                        t_total=0.0_dp, &
-                                        t_write=0.0_dp
-       CHARACTER(LEN=*), PARAMETER  ::  fmt='(F6.2)'
+       CHARACTER(LEN=*), INTENT(IN) ::  mode   !< name of the calling mode
+       CHARACTER(LEN=*), INTENT(IN) ::  budget !< name of the timing budget
+
+       REAL(dp), SAVE ::  t0               !< begin of timing interval
+       REAL(dp), SAVE ::  t1               !< end of timing interval
+       REAL(dp), SAVE ::  t_comp  = 0.0_dp !< computation timing budget
+       REAL(dp), SAVE ::  t_alloc = 0.0_dp !< allocation timing budget
+       REAL(dp), SAVE ::  t_init  = 0.0_dp !< initialization timing budget
+       REAL(dp), SAVE ::  t_read  = 0.0_dp !< reading timing budget
+       REAL(dp), SAVE ::  t_total = 0.0_dp !< total time
+       REAL(dp), SAVE ::  t_write = 0.0_dp !< writing timing budget
+
+       CHARACTER(LEN=*), PARAMETER  ::  fmt='(F6.2)' !< floating-point output format
 
 
Index: palm/trunk/UTIL/inifor/src/inifor_defs.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_defs.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor_defs.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+! 
+! 
+! 3537 2018-11-20 10:53:14Z eckhard
 ! Bumped version number
 ! 
@@ -59,5 +63,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -68,6 +72,7 @@
  
  IMPLICIT NONE
- 
- ! Parameters for type definitions
+
+! 
+!-- Parameters for type definitions
  INTEGER, PARAMETER  ::  dp    = 8   !< double precision (8 bytes = 64 bits)
  INTEGER, PARAMETER  ::  sp    = 4   !< single precision (4 bytes = 32 bits)
@@ -78,29 +83,38 @@
  INTEGER, PARAMETER  ::  DATE  = 10  !< length of date strings
 
- ! Trigonomentry
+! 
+!-- Trigonomentry
  REAL(dp), PARAMETER ::  PI = 3.14159265358979323846264338_dp !< Ratio of a circle's circumference to its diamter [-]
  REAL(dp), PARAMETER ::  TO_RADIANS = PI / 180.0_dp           !< Conversion factor from degrees to radiant [-]
  REAL(dp), PARAMETER ::  TO_DEGREES = 180.0_dp / PI           !< Conversion factor from radians to degrees [-]
 
- ! COSMO-DE parameters
- INTEGER, PARAMETER  ::  WATER_ID = 9                         !< Integer corresponding to the water soil type in COSMO-DE [-]
- REAL(dp), PARAMETER ::  EARTH_RADIUS = 6371229.0_dp          !< Earth radius used in COSMO-DE [m]
- REAL(dp), PARAMETER ::  P_SL = 1e5_dp                        !< Reference pressure for computation of COSMO-DE's basic state pressure [Pa]
- REAL(dp), PARAMETER ::  T_SL = 288.15_dp                     !< Reference temperature for computation of COSMO-DE's basic state pressure [K]
- REAL(dp), PARAMETER ::  BETA = 42.0_dp                       !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic state pressure [K]
- REAL(dp), PARAMETER ::  RD   = 287.05_dp                     !< specific gas constant of dry air, used in computation of COSMO-DE's basic state [J/kg/K]
- REAL(dp), PARAMETER ::  RV   = 461.51_dp                     !< specific gas constant of water vapor [J/kg/K]
- REAL(dp), PARAMETER ::  G    = 9.80665_dp                    !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s]
- REAL(dp), PARAMETER ::  RHO_L = 1e3_dp                       !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3], in [kg/m^3]
- REAL(dp), PARAMETER ::  HECTO = 100_dp                       !< unit conversion factor from hPa to Pa
+! 
+!-- COSMO-DE parameters
+ INTEGER, PARAMETER  ::  WATER_ID = 9                !< Integer corresponding to the water soil type in COSMO-DE [-]
+ REAL(dp), PARAMETER ::  EARTH_RADIUS = 6371229.0_dp !< Earth radius used in COSMO-DE [m]
+ REAL(dp), PARAMETER ::  P_SL = 1e5_dp               !< Reference pressure for computation of COSMO-DE's basic state pressure [Pa]
+ REAL(dp), PARAMETER ::  T_SL = 288.15_dp            !< Reference temperature for computation of COSMO-DE's basic state pressure [K]
+ REAL(dp), PARAMETER ::  BETA = 42.0_dp              !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic
+                                                     !< state pressure [K]
+ REAL(dp), PARAMETER ::  RD   = 287.05_dp            !< specific gas constant of dry air, used in computation of COSMO-DE's basic
+                                                     !< state [J/kg/K]
+ REAL(dp), PARAMETER ::  RV   = 461.51_dp            !< specific gas constant of water vapor [J/kg/K]
+ REAL(dp), PARAMETER ::  G    = 9.80665_dp           !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic
+                                                     !< state [m/s/s]
+ REAL(dp), PARAMETER ::  RHO_L = 1e3_dp              !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3],
+                                                     !< in [kg/m^3]
+ 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]
- REAL(dp), PARAMETER ::  CP_PALM = 1005.0_dp                  !< heat capacity of dry air at constant pressure, used in computation of PALM-4U's potential temperature [J/kg/K]
+!
+!-- 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]
+ REAL(dp), PARAMETER ::  CP_PALM = 1005.0_dp  !< heat capacity of dry air at constant pressure, used in computation of PALM-4U's potential temperature [J/kg/K]
 
- ! INIFOR parameters
- INTEGER, PARAMETER          ::  FILL_ITERATIONS = 5          !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-]
+!
+!-- INIFOR parameters
+ INTEGER, PARAMETER          ::  FILL_ITERATIONS = 5          !< Number of iterations for extrapolating soil data into COSMO-DE
+                                                              !< water cells [-]
  INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
  REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
Index: palm/trunk/UTIL/inifor/src/inifor_grid.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_grid.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor_grid.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+!
+! 
+! 3537 2018-11-20 10:53:14Z eckhard
 ! Read COSMO domain extents and soil depths from input files
 ! Report averaging mode and debugging mode in log
@@ -69,5 +73,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 !------------------------------------------------------------------------------!
@@ -193,101 +197,109 @@
     INTEGER ::  step_hour !< number of hours between forcing time steps
 
-    LOGICAL ::  init_variables_required
-    LOGICAL ::  boundary_variables_required
-    LOGICAL ::  ls_forcing_variables_required
-    LOGICAL ::  profile_grids_required
-    LOGICAL ::  surface_forcing_required
+    LOGICAL ::  init_variables_required       !< flag controlling whether init variables are to be processed
+    LOGICAL ::  boundary_variables_required   !< flag controlling whether boundary grids are to be allocated and boundary variables are to be computed
+    LOGICAL ::  ls_forcing_variables_required !< flag controlling whether large-scale forcing variables are to be computed
+    LOGICAL ::  surface_forcing_required      !< flag controlling whether surface forcing variables are to be computed
 
     TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
     TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
-    TYPE(nc_var)                      ::  cosmo_var           !< COSMO-DE dummy variable, used for reading HHL, rlon, rlat
+    TYPE(nc_var)                      ::  cosmo_var           !< COSMO dummy variable, used for reading HHL, rlon, rlat
 
     TYPE(grid_definition), TARGET ::  palm_grid                       !< PALM-4U grid in the target system (COSMO-DE rotated-pole)
     TYPE(grid_definition), TARGET ::  palm_intermediate               !< PALM-4U grid with coarse vertical grid wiht levels interpolated from COSMO-DE grid
     TYPE(grid_definition), TARGET ::  cosmo_grid                      !< target system (COSMO-DE rotated-pole)
-    TYPE(grid_definition), TARGET ::  scalars_east_grid               !< 
-    TYPE(grid_definition), TARGET ::  scalars_west_grid               !< 
-    TYPE(grid_definition), TARGET ::  scalars_north_grid              !< 
-    TYPE(grid_definition), TARGET ::  scalars_south_grid              !< 
-    TYPE(grid_definition), TARGET ::  scalars_top_grid                !< 
-    TYPE(grid_definition), TARGET ::  scalars_east_intermediate       !< 
-    TYPE(grid_definition), TARGET ::  scalars_west_intermediate       !< 
-    TYPE(grid_definition), TARGET ::  scalars_north_intermediate      !< 
-    TYPE(grid_definition), TARGET ::  scalars_south_intermediate      !< 
-    TYPE(grid_definition), TARGET ::  scalars_top_intermediate        !< 
-    TYPE(grid_definition), TARGET ::  u_initial_grid                  !< 
-    TYPE(grid_definition), TARGET ::  u_east_grid                     !< 
-    TYPE(grid_definition), TARGET ::  u_west_grid                     !< 
-    TYPE(grid_definition), TARGET ::  u_north_grid                    !< 
-    TYPE(grid_definition), TARGET ::  u_south_grid                    !< 
-    TYPE(grid_definition), TARGET ::  u_top_grid                      !< 
-    TYPE(grid_definition), TARGET ::  u_initial_intermediate          !<
-    TYPE(grid_definition), TARGET ::  u_east_intermediate             !< 
-    TYPE(grid_definition), TARGET ::  u_west_intermediate             !< 
-    TYPE(grid_definition), TARGET ::  u_north_intermediate            !< 
-    TYPE(grid_definition), TARGET ::  u_south_intermediate            !< 
-    TYPE(grid_definition), TARGET ::  u_top_intermediate              !< 
-    TYPE(grid_definition), TARGET ::  v_initial_grid                  !< 
-    TYPE(grid_definition), TARGET ::  v_east_grid                     !< 
-    TYPE(grid_definition), TARGET ::  v_west_grid                     !< 
-    TYPE(grid_definition), TARGET ::  v_north_grid                    !< 
-    TYPE(grid_definition), TARGET ::  v_south_grid                    !< 
-    TYPE(grid_definition), TARGET ::  v_top_grid                      !< 
-    TYPE(grid_definition), TARGET ::  v_initial_intermediate          !<
-    TYPE(grid_definition), TARGET ::  v_east_intermediate             !< 
-    TYPE(grid_definition), TARGET ::  v_west_intermediate             !< 
-    TYPE(grid_definition), TARGET ::  v_north_intermediate            !< 
-    TYPE(grid_definition), TARGET ::  v_south_intermediate            !< 
-    TYPE(grid_definition), TARGET ::  v_top_intermediate              !< 
-    TYPE(grid_definition), TARGET ::  w_initial_grid                  !< 
-    TYPE(grid_definition), TARGET ::  w_east_grid                     !< 
-    TYPE(grid_definition), TARGET ::  w_west_grid                     !< 
-    TYPE(grid_definition), TARGET ::  w_north_grid                    !< 
-    TYPE(grid_definition), TARGET ::  w_south_grid                    !< 
-    TYPE(grid_definition), TARGET ::  w_top_grid                      !< 
-    TYPE(grid_definition), TARGET ::  w_initial_intermediate          !<
-    TYPE(grid_definition), TARGET ::  w_east_intermediate             !< 
-    TYPE(grid_definition), TARGET ::  w_west_intermediate             !< 
-    TYPE(grid_definition), TARGET ::  w_north_intermediate            !< 
-    TYPE(grid_definition), TARGET ::  w_south_intermediate            !< 
-    TYPE(grid_definition), TARGET ::  w_top_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(grid_definition), TARGET ::  scalars_east_grid               !< grid for eastern scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_west_grid               !< grid for western scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_north_grid              !< grid for northern scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_south_grid              !< grid for southern scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_top_grid                !< grid for top scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_east_intermediate       !< intermediate grid for eastern scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_west_intermediate       !< intermediate grid for western scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_north_intermediate      !< intermediate grid for northern scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_south_intermediate      !< intermediate grid for southern scalar boundary condition
+    TYPE(grid_definition), TARGET ::  scalars_top_intermediate        !< intermediate grid for top scalar boundary condition
+    TYPE(grid_definition), TARGET ::  u_initial_grid                  !< grid for u initial condition
+    TYPE(grid_definition), TARGET ::  u_east_grid                     !< grid for eastern u boundary condition
+    TYPE(grid_definition), TARGET ::  u_west_grid                     !< grid for western u boundary condition
+    TYPE(grid_definition), TARGET ::  u_north_grid                    !< grid for northern u boundary condition
+    TYPE(grid_definition), TARGET ::  u_south_grid                    !< grid for southern u boundary condition
+    TYPE(grid_definition), TARGET ::  u_top_grid                      !< grid for top u boundary condition
+    TYPE(grid_definition), TARGET ::  u_initial_intermediate          !< intermediate grid for u initial condition
+    TYPE(grid_definition), TARGET ::  u_east_intermediate             !< intermediate grid for eastern u boundary condition
+    TYPE(grid_definition), TARGET ::  u_west_intermediate             !< intermediate grid for western u boundary condition
+    TYPE(grid_definition), TARGET ::  u_north_intermediate            !< intermediate grid for northern u boundary condition
+    TYPE(grid_definition), TARGET ::  u_south_intermediate            !< intermediate grid for southern u boundary condition
+    TYPE(grid_definition), TARGET ::  u_top_intermediate              !< intermediate grid for top u boundary condition
+    TYPE(grid_definition), TARGET ::  v_initial_grid                  !< grid for v initial condition
+    TYPE(grid_definition), TARGET ::  v_east_grid                     !< grid for eastern v boundary condition
+    TYPE(grid_definition), TARGET ::  v_west_grid                     !< grid for western v boundary condition
+    TYPE(grid_definition), TARGET ::  v_north_grid                    !< grid for northern v boundary condition
+    TYPE(grid_definition), TARGET ::  v_south_grid                    !< grid for southern v boundary condition
+    TYPE(grid_definition), TARGET ::  v_top_grid                      !< grid for top v boundary condition
+    TYPE(grid_definition), TARGET ::  v_initial_intermediate          !< intermediate grid for v initial condition
+    TYPE(grid_definition), TARGET ::  v_east_intermediate             !< intermediate grid for eastern v boundary condition
+    TYPE(grid_definition), TARGET ::  v_west_intermediate             !< intermediate grid for western v boundary condition
+    TYPE(grid_definition), TARGET ::  v_north_intermediate            !< intermediate grid for northern v boundary condition
+    TYPE(grid_definition), TARGET ::  v_south_intermediate            !< intermediate grid for southern v boundary condition
+    TYPE(grid_definition), TARGET ::  v_top_intermediate              !< intermediate grid for top v boundary condition
+    TYPE(grid_definition), TARGET ::  w_initial_grid                  !< grid for w initial condition
+    TYPE(grid_definition), TARGET ::  w_east_grid                     !< grid for eastern w boundary condition
+    TYPE(grid_definition), TARGET ::  w_west_grid                     !< grid for western w boundary condition
+    TYPE(grid_definition), TARGET ::  w_north_grid                    !< grid for northern w boundary condition
+    TYPE(grid_definition), TARGET ::  w_south_grid                    !< grid for southern w boundary condition
+    TYPE(grid_definition), TARGET ::  w_top_grid                      !< grid for top w boundary condition
+    TYPE(grid_definition), TARGET ::  w_initial_intermediate          !< intermediate grid for w initial condition
+    TYPE(grid_definition), TARGET ::  w_east_intermediate             !< intermediate grid for eastern w boundary condition
+    TYPE(grid_definition), TARGET ::  w_west_intermediate             !< intermediate grid for western w boundary condition
+    TYPE(grid_definition), TARGET ::  w_north_intermediate            !< intermediate grid for northern w boundary condition
+    TYPE(grid_definition), TARGET ::  w_south_intermediate            !< intermediate grid for southern w boundary condition
+    TYPE(grid_definition), TARGET ::  w_top_intermediate              !< intermediate grid for top w boundary condition
+    TYPE(grid_definition), TARGET ::  north_averaged_scalar_profile   !< grid of the northern scalar averaging region
+    TYPE(grid_definition), TARGET ::  south_averaged_scalar_profile   !< grid of the southern scalar averaging region
+    TYPE(grid_definition), TARGET ::  west_averaged_scalar_profile    !< grid of the western scalar averaging region
+    TYPE(grid_definition), TARGET ::  east_averaged_scalar_profile    !< grid of the eastern scalar averaging region
+    TYPE(grid_definition), TARGET ::  averaged_scalar_profile         !< grid of the central scalar averaging region
+    TYPE(grid_definition), TARGET ::  averaged_w_profile              !< grid of the central w-velocity averaging region
 
     TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
  
     NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude,             &
-                      dz_max, dz_stretch_factor, dz_stretch_level,             & !< old grid stretching parameters
-                      dz_stretch_level_start, dz_stretch_level_end               !< new grid stretching parameters
+                      dz_max, dz_stretch_factor, dz_stretch_level,             &
+                      dz_stretch_level_start, dz_stretch_level_end
     NAMELIST /d3par/  end_time
     
-    CHARACTER(LEN=LNAME) ::  nc_source_text     = ''  !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...'
-    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  flow_files
-    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_moisture_files
-    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_files
-    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'
-    CHARACTER(LEN=SNAME) ::  radiation_suffix     !< Suffix of radiation input files, e.g. 'radiation'
-    CHARACTER(LEN=SNAME) ::  soilmoisture_suffix  !< Suffix of input files for soil moisture spin-up, e.g. 'soilmoisture'
+    CHARACTER(LEN=LNAME) ::  nc_source_text = ''  !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...'
+
+    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  flow_files          !< list of atmospheric input files (<prefix>YYYYMMDDHH-flow.nc)
+    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_moisture_files !< list of precipitation input files (<prefix>YYYYMMDDHH-soilmoisture.nc) 
+    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_files          !< list of soil input files (temperature, moisture, <prefix>YYYYMMDDHH-soil.nc)
+    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  radiation_files     !< list of radiation input files (<prefix>YYYYMMDDHH-rad.nc)
+
+    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'
+    CHARACTER(LEN=SNAME) ::  radiation_suffix     !< suffix of radiation input files, e.g. 'radiation'
+    CHARACTER(LEN=SNAME) ::  soilmoisture_suffix  !< suffix of input files for soil moisture spin-up, e.g. 'soilmoisture'
                           
-    TYPE(nc_file) ::  output_file
-
-    TYPE(inifor_config) ::  cfg !< Container of INIFOR configuration
+    TYPE(nc_file) ::  output_file !< metadata of the dynamic driver 
+
+    TYPE(inifor_config) ::  cfg !< container of the INIFOR command-line configuration
 
  CONTAINS
     
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine initializes INIFOR. This includes parsing command-line options,
+!> setting the names of the input and output files, reading INIFOR's namelist
+!> file as well as reading in and setting grid parameters for defining
+!> interplation grids later in setup_grids().
+!------------------------------------------------------------------------------!
     SUBROUTINE setup_parameters()
-       INTEGER ::  k
+       INTEGER ::  k !< loop index
 
 !
@@ -300,20 +312,24 @@
        start_hour_soilmoisture = - (4 * 7 * 24) - 2
 
-       ! COSMO-DE default rotated pole
+!
+!--    COSMO-DE and -D2 rotated pole position
        phi_n     =   40.0_dp * TO_RADIANS
        phi_equat =   50.0_dp * TO_RADIANS
        lambda_n  = -170.0_dp * TO_RADIANS
 
-       ! Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
+!
+!--    Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
        origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
        origin_lon = 13.082744_dp * TO_RADIANS
        cfg % z0 = 35.0_dp
 
-       ! Default atmospheric parameters
+!
+!--    Default atmospheric parameters
        cfg % ug = 0.0_dp
        cfg % vg = 0.0_dp
        cfg % p0 = P_SL
 
-       ! Parameters for file names
+!
+!--    Parameters for file names
        start_hour_flow = 0
        start_hour_soil = 0
@@ -323,5 +339,5 @@
        step_hour = FORCING_STEP
 
-       input_prefix = 'laf'  ! 'laf' for COSMO-DE analyses
+       input_prefix = 'laf'
        cfg % flow_prefix = input_prefix
        cfg % input_prefix = input_prefix
@@ -342,5 +358,6 @@
 !------------------------------------------------------------------------------
 
-       ! Set default paths and modes
+!
+!--    Set default paths and modes
        cfg % input_path         = './'
        cfg % hhl_file           = ''
@@ -353,5 +370,6 @@
        cfg % averaging_mode = 'level'
 
-       ! Overwrite defaults with user configuration
+!
+!--    Overwrite defaults with user configuration
        CALL parse_command_line_arguments( cfg )
 
@@ -380,6 +398,6 @@
        END IF
 
-
-       ! Set default file paths, if not specified by user.
+!
+!--    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'
@@ -399,5 +417,6 @@
 
  CALL run_control('time', 'init')
-       ! Read in namelist parameters
+ !
+ !--   Read in namelist parameters
        OPEN(10, FILE=cfg % namelist_file)
        READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz
@@ -408,5 +427,6 @@
        end_hour = CEILING( end_time / 3600.0 * step_hour )
 
-       ! Generate input file lists
+!
+!--    Generate input file lists
        CALL get_input_file_list(                                               &
           cfg % start_date, start_hour_flow, end_hour, step_hour,              &
@@ -437,5 +457,6 @@
 
  CALL run_control('time', 'init')
-       ! Read COSMO-DE soil type map
+!
+!--    Read COSMO soil type map
        cosmo_var % name = 'SOILTYP'
        CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp)
@@ -465,6 +486,6 @@
        CALL report('setup_parameters', message)
 
-
-       ! Read in COSMO-DE heights of half layers (vertical cell faces)
+!
+!--    Read in COSMO heights of half layers (vertical cell faces)
        cosmo_var % name = 'HHL'
        CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl)
@@ -486,5 +507,6 @@
  CALL run_control('time', 'comp')
 
-       ! Appoximate COSMO-DE heights of full layers (cell centres)
+!
+!--    Appoximate COSMO-DE heights of full layers (cell centres)
        ALLOCATE( hfl(nlon, nlat, nlev-1) )
        ALLOCATE( d_depth(ndepths), d_depth_rho_inv(ndepths) )
@@ -494,5 +516,6 @@
        d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L )
 
-       ! Appoximate COSMO-DE heights of full layers (cell centres)
+!
+!--    Appoximate COSMO-DE heights of full layers (cell centres)
        DO k = 1, nlev-1
           hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) +                                 &
@@ -506,13 +529,16 @@
 ! Section 4.2: PALM-4U parameters
 !------------------------------------------------------------------------------
-       ! PALM-4U domain extents
+!
+!--    PALM-4U domain extents
        lx = (nx+1) * dx
        ly = (ny+1) * dy
        
-       ! PALM-4U point of Earth tangency
+!
+!--    PALM-4U point of Earth tangency
        x0 = 0.0_dp
        y0 = 0.0_dp
 
-       ! time vector
+!
+!--    time vector
        nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1
        ALLOCATE( time(nt) )
@@ -521,5 +547,6 @@
  CALL run_control('time', 'init')
 
-! Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid
+!
+!--    Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid
        phi_c    = TO_RADIANS *                                                 &
                   phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
@@ -530,13 +557,15 @@
                               0.0_dp )
 
-! Set gamma according to whether PALM domain is in the northern or southern
-! hemisphere of the COSMO-DE rotated-pole system. Gamma assumes either the
-! value 0 or PI and is needed to work around around a bug in the rotated-pole
-! coordinate transformations.
+!
+!--    Set gamma according to whether PALM domain is in the northern or southern
+!--    hemisphere of the COSMO rotated-pole system. Gamma assumes either the
+!--    value 0 or PI and is needed to work around around a bug in the
+!--    rotated-pole coordinate transformations.
        gam = gamma_from_hemisphere(origin_lat, phi_equat)
 
-! Compute the north pole of the rotated-pole grid centred at the PALM-4U domain
-! centre. The resulting (phi_cn, lambda_cn) are coordinates in COSMO-DE's
-! rotated-pole grid.
+!
+!--    Compute the north pole of the rotated-pole grid centred at the PALM-4U
+!--    domain centre. The resulting (phi_cn, lambda_cn) are coordinates in
+!--    COSMO-DE's rotated-pole grid.
        phi_cn    = phic_to_phin(phi_c) 
        lambda_cn = lamc_to_lamn(phi_c, lambda_c) 
@@ -573,5 +602,6 @@
 !------------------------------------------------------------------------------
 
-    ! Compute coordiantes of the PALM centre in the source (COSMO) system
+!
+!-- 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,             &
@@ -595,6 +625,6 @@
     CALL report( 'setup_parameters', message )
 
-
-    ! Compute boundaries of the central averaging box
+!
+!-- Compute boundaries of the central averaging box
     averaging_angle = cfg % averaging_angle * TO_RADIANS
     lam_east = lam_centre + 0.5_dp * averaging_angle
@@ -605,5 +635,6 @@
     averaging_width_ns = averaging_angle * EARTH_RADIUS
 
-    ! Coriolis parameter
+!
+!-- Coriolis parameter
     f3 = 2.0_dp * OMEGA * SIN(                                                 &
        TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,&
@@ -618,7 +649,8 @@
 ! Description:
 ! ------------
-!> Initializes the COSMO-DE, PALM-4U, PALM-4U boundary grids.
+!> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in partucular their
+!> coordinates and interpolation weights
 !------------------------------------------------------------------------------!
-    SUBROUTINE setup_grids() ! setup_grids(inifor_settings(with nx, ny, nz,...))
+    SUBROUTINE setup_grids()
        CHARACTER ::  interp_mode
 
@@ -626,8 +658,8 @@
 ! Section 0: Define base PALM-4U grid coordinate vectors
 !------------------------------------------------------------------------------
-       ! palm x y z, we allocate the column to nz+1 in order to include the top
-       ! scalar boundary. The interpolation grids will be associated with
-       ! a shorter column that omits the top element.
-       
+!
+!--    palm x y z, we allocate the column to nz+1 in order to include the top
+!--    scalar boundary. The interpolation grids will be associated with
+!--    a shorter column that omits the top element.
        ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
        CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x)
@@ -642,6 +674,7 @@
        z_top = z_column(nz+1)
 
-       ! palm xu yv zw, compared to the scalar grid, velocity coordinates
-       ! contain one element less.
+!
+!--    palm xu yv zw, compared to the scalar grid, velocity coordinates
+!--    contain one element less.
        ALLOCATE( xu(1:nx),  yv(1:ny), zw(1:nz-1), zw_column(1:nz))
        CALL linspace(dx, lx - dx, xu)
@@ -661,5 +694,6 @@
                nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode)
 
-       ! Subtracting 1 because arrays will be allocated with nlon + 1 elements.
+!
+!--    Subtracting 1 because arrays will be allocated with nlon + 1 elements.
        CALL init_grid_definition('cosmo-de', grid=cosmo_grid,                  &
                xmin=lonmin, xmax=lonmax,                                       &
@@ -668,8 +702,9 @@
                nx=nlon-1, ny=nlat-1, nz=nlev-1)
 
-! Define intermediate grid. This is the same as palm_grid except with a much
-! coarser vertical grid. The vertical levels are interpolated in each PALM-4U
-! column from COSMO-DE's secondary levels. The main levels are then computed as
-! the averages of the bounding secondary levels.
+!
+!--    Define intermediate grid. This is the same as palm_grid except with a
+!--    much coarser vertical grid. The vertical levels are interpolated in each
+!--    PALM column from COSMO's secondary levels. The main levels are then
+!--    computed as the averages of the bounding secondary levels.
        CALL init_grid_definition('palm intermediate', grid=palm_intermediate,  &
                xmin=0.0_dp, xmax=lx,                                           &
@@ -1007,35 +1042,4 @@
                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
-       !   ! 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
 
 !                                                                              
@@ -1091,4 +1095,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Driver for computing neighbour indices and weights for horizontal and
+!> vertical interpolation.
+!------------------------------------------------------------------------------!
     SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
 
@@ -1106,20 +1116,25 @@
 ! Section 1: Horizontal interpolation                                        
 !------------------------------------------------------------------------------
-       ! Select horizontal coordinates according to kind of points (s/w, u, v)
+!
+!--    Select horizontal coordinates according to kind of points (s/w, u, v)
        SELECT CASE(kind)
 
-       CASE('s') ! scalars
+!
+!--    scalars
+       CASE('s')
 
           cosmo_lat => cosmo_grid % lat
           cosmo_lon => cosmo_grid % lon
           cosmo_h   => cosmo_grid % hfl
-
-       CASE('w') ! vertical velocity
+!
+!--    vertical velocity
+       CASE('w')
 
           cosmo_lat => cosmo_grid % lat
           cosmo_lon => cosmo_grid % lon
           cosmo_h   => cosmo_grid % hhl
-
-       CASE('u') ! x velocity
+!
+!--    x velocity
+       CASE('u')
 
           cosmo_lat => cosmo_grid % lat
@@ -1127,5 +1142,7 @@
           cosmo_h   => cosmo_grid % hfl
 
-       CASE('v') ! y velocity
+!
+!--    y velocity
+       CASE('v')
 
           cosmo_lat => cosmo_grid % latv
@@ -1153,10 +1170,11 @@
 !------------------------------------------------------------------------------
 
-       ! 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 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.
@@ -1169,6 +1187,7 @@
           intermediate_grid % h(:,:,:) = - EARTH_RADIUS
 
-          ! For w points, use hhl, for scalars use hfl
-          ! compute the full heights for the intermediate grids
+!
+!--       For w points, use hhl, for scalars use hfl
+!--       compute the full heights for the intermediate grids
           CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid)
           CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
@@ -1246,5 +1265,6 @@
            grid % z => z
 
-           ! Allocate neighbour indices and weights
+!
+!--        Allocate neighbour indices and weights
            IF (TRIM(ic_mode) .NE. 'profile')  THEN
               ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
@@ -1274,5 +1294,6 @@
            )
 
-           ! Allocate neighbour indices and weights
+!
+!--        Allocate neighbour indices and weights
            ALLOCATE( grid % ii(0:nx, 0:ny, 4),                                 &
                      grid % jj(0:nx, 0:ny, 4) )
@@ -1283,6 +1304,7 @@
            grid % w_horiz(:,:,:)   = 0.0_dp
         
-        ! This mode initializes a Cartesian PALM-4U grid and adds the
-        ! corresponding latitudes and longitudes of the rotated pole grid.
+!
+!--     This mode initializes a Cartesian PALM-4U grid and adds the
+!--     corresponding latitudes and longitudes of the rotated pole grid.
         CASE('palm')
 
@@ -1301,6 +1323,7 @@
            grid % name(3) = 'z'
 
-           !TODO: Remove use of global dx, dy, dz variables. Consider
-           !TODO:   associating global x,y, and z arrays.
+!
+!--        TODO: Remove use of global dx, dy, dz variables. Consider
+!--        TODO: associating global x,y, and z arrays.
            ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
            ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
@@ -1314,5 +1337,6 @@
            grid % depths => depths
 
-           ! Allocate neighbour indices and weights
+!
+!--        Allocate neighbour indices and weights
            IF (TRIM(ic_mode) .NE. 'profile')  THEN
               ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
@@ -1329,6 +1353,7 @@
            grid % name(3) = 'interpolated hhl or hfl'
 
-           !TODO: Remove use of global dx, dy, dz variables. Consider
-           !TODO:   associating global x,y, and z arrays.
+!
+!--        TODO: Remove use of global dx, dy, dz variables. Consider
+!--        TODO: associating global x,y, and z arrays.
            ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
            ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
@@ -1340,11 +1365,13 @@
            grid % depths => depths
 
-           ! Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude
+!
+!--        Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude
            ALLOCATE( grid % clon(0:nx, 0:ny),   grid % clat(0:nx, 0:ny)  )
            ALLOCATE( grid % clonu(1:nx, 0:ny),  grid % clatu(1:nx, 0:ny) )
            ALLOCATE( grid % clonv(0:nx, 1:ny),  grid % clatv(0:nx, 1:ny) )
 
-           ! Compute rotated-pole coordinates of...
-           ! ... PALM-4U centres
+!
+!--        Compute rotated-pole coordinates of...
+!--        ... PALM-4U centres
            CALL rotate_to_cosmo(                                               &
               phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
@@ -1356,5 +1383,6 @@
            )
 
-           ! ... PALM-4U u winds
+!
+!--        ... PALM-4U u winds
            CALL rotate_to_cosmo(                                               &
               phir = project( grid % y,  y0, EARTH_RADIUS ), & ! = plate-carree latitude
@@ -1366,5 +1394,6 @@
            )
 
-           ! ... PALM-4U v winds
+!
+!--        ... PALM-4U v winds
            CALL rotate_to_cosmo(                                               &
               phir = project( grid % yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude
@@ -1376,5 +1405,6 @@
            )
 
-           ! Allocate neighbour indices and weights
+!
+!--        Allocate neighbour indices and weights
            ALLOCATE( grid % ii(0:nx, 0:ny, 4),                                 &
                      grid % jj(0:nx, 0:ny, 4) )
@@ -1398,6 +1428,7 @@
            grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny)
 
-           ! Point to heights of half levels (hhl) and compute heights of full
-           ! levels (hfl) as arithmetic averages
+!
+!--        Point to heights of half levels (hhl) and compute heights of full
+!--        levels (hfl) as arithmetic averages
            grid % hhl => hhl
            grid % hfl => hfl
@@ -1472,11 +1503,13 @@
        avg_grid % kind = TRIM(kind)
 
-       ! Find and store COSMO columns that fall into the coordinate range
-       ! given by avg_grid % clon, % clat
+!
+!--    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
+!
+!--    Compute average COSMO levels in the averaging region
        SELECT CASE(avg_grid % kind)
 
@@ -1494,5 +1527,6 @@
        END SELECT
 
-       ! For level-besed averaging, compute average heights
+!
+!--    For level-besed averaging, compute average heights
        level_based_averaging = .TRUE.
        IF (level_based_averaging)  THEN
@@ -1503,7 +1537,8 @@
        END IF
 
-       ! Copy averaged grid to all COSMO columns, leads to computing the same
-       ! vertical interpolation weights for all columns and to COSMO grid level
-       ! based averaging onto the averaged COSMO heights.
+!
+!--    Copy averaged grid to all COSMO columns, 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(                                                              &
@@ -1513,5 +1548,6 @@
        END IF
 
-       ! Compute vertical weights and neighbours
+!
+!--    Compute vertical weights and neighbours
        CALL find_vertical_neighbours_and_weights_average(avg_grid)
 
@@ -1550,6 +1586,6 @@
        END SELECT
 
-
-       ! FIXME: make dlon, dlat parameters of the grid_defintion type
+!
+!--    FIXME: make dlon, dlat parameters of the grid_defintion type
        dlon = cosmo_lon(1) - cosmo_lon(0)
        dlat = cosmo_lat(1) - cosmo_lat(0)
@@ -1610,8 +1646,9 @@
        REAL(dp) ::  dz_level_end, dz_stretched
 
-       INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
-       INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
+       INTEGER ::  dz_stretch_level_end_index(9)      !< vertical grid level index until which the vertical grid spacing is stretched
+       INTEGER ::  dz_stretch_level_start_index(9)    !< vertical grid level index above which the vertical grid spacing is stretched
        INTEGER ::  dz_stretch_level_index = 0
        INTEGER ::  k, n, number_dz
+
 !
 !-- Compute height of u-levels from constant grid length and dz stretch factors
@@ -1826,4 +1863,5 @@
 
 
+!------------------------------------------------------------------------------!
 ! Description: [PALM subroutine]
 ! -----------------------------------------------------------------------------!
@@ -1936,8 +1974,8 @@
              
 !
-!--                stretch_factor_1 is taken to guarantee that the stretching
-!--                procedure ends as close as possible to dz_stretch_level_end.
-!--                stretch_factor_2 would guarantee that the stretched dz(n) is
-!--                equal to dz(n+1) after l_rounded grid levels.
+!--          stretch_factor_1 is taken to guarantee that the stretching
+!--          procedure ends as close as possible to dz_stretch_level_end.
+!--          stretch_factor_2 would guarantee that the stretched dz(n) is
+!--          equal to dz(n+1) after l_rounded grid levels.
              IF (delta_total_new < delta_total_old) THEN
                 dz_stretch_factor_array(n) = stretch_factor_1
@@ -1974,4 +2012,10 @@
  END SUBROUTINE calculate_stretching_factor
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computes the midpoint between two neighbouring coordinates in the given
+!> coordinate vector 'z' and stores it in 'zw'.
+!------------------------------------------------------------------------------!
     SUBROUTINE midpoints(z, zw)
 
@@ -1987,5 +2031,9 @@
     END SUBROUTINE midpoints
 
-
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Defines INFOR's IO groups.
+!------------------------------------------------------------------------------!
     SUBROUTINE setup_io_groups()
 
@@ -1994,6 +2042,7 @@
        ngroups = 16
        ALLOCATE( io_group_list(ngroups) )
-       
-       !soil temp
+
+!
+!--    soil temp
        io_group_list(1) = init_io_group(                                       &
           in_files = soil_files,                                               &
@@ -2003,5 +2052,6 @@
        )
 
-       !soil water
+!
+!--    soil water
        io_group_list(2) = init_io_group(                                       &
           in_files = soil_files,                                               &
@@ -2011,6 +2061,7 @@
        )
 
-       !potential temperature, surface pressure, specific humidity including
-       !nudging and subsidence, and geostrophic winds ug, vg
+!
+!--    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,                                               &
@@ -2025,5 +2076,6 @@
        )
 
-       !Moved to therodynamic io_group
+!
+!--    Moved to therodynamic io_group
        !io_group_list(4) = init_io_group(                                       &
        !   in_files = flow_files,                                               &
@@ -2033,5 +2085,6 @@
        !)
 
-       !u and v velocity, ug anv vg moved to thermodynamic io group. 
+!
+!--    u and v velocity
        io_group_list(5) = init_io_group(                                       &
           in_files = flow_files,                                               &
@@ -2043,5 +2096,6 @@
        )
    
-       !!v velocity, deprecated!
+!
+!--    v velocity, deprecated!
        !io_group_list(6) = init_io_group(                                       &
        !   in_files = flow_files,                                               &
@@ -2052,5 +2106,6 @@
        !io_group_list(6) % to_be_processed = .FALSE.
    
-       !w velocity and subsidence and w nudging
+!
+!--    w velocity and subsidence and w nudging
        io_group_list(7) = init_io_group(                                       &
           in_files = flow_files,                                               &
@@ -2059,6 +2114,6 @@
           kind = 'scalar'                                                      &
        )
-
-       !rain
+!
+!--    rain
        io_group_list(8) = init_io_group(                                       &
           in_files = soil_moisture_files,                                      &
@@ -2068,6 +2123,6 @@
        )
        io_group_list(8) % to_be_processed = .FALSE.
-
-       !snow
+!
+!--    snow
        io_group_list(9) = init_io_group(                                       &
           in_files = soil_moisture_files,                                      &
@@ -2077,6 +2132,6 @@
        )
        io_group_list(9) % to_be_processed = .FALSE.
-
-       !graupel
+!
+!--    graupel
        io_group_list(10) = init_io_group(                                      &
           in_files = soil_moisture_files,                                      &
@@ -2086,6 +2141,6 @@
        )
        io_group_list(10) % to_be_processed = .FALSE.
-
-       !evapotranspiration
+!
+!--    evapotranspiration
        io_group_list(11) = init_io_group(                                      &
           in_files = soil_moisture_files,                                      &
@@ -2095,6 +2150,6 @@
        )
        io_group_list(11) % to_be_processed = .FALSE.
-
-       !2m air temperature
+!
+!--    2m air temperature
        io_group_list(12) = init_io_group(                                      &
           in_files = soil_moisture_files,                                      &
@@ -2103,6 +2158,6 @@
           kind = 'surface'                                                     &
        )
-
-       !incoming diffusive sw flux
+!
+!--    incoming diffusive sw flux
        io_group_list(13) = init_io_group(                                      &
           in_files = radiation_files,                                          &
@@ -2112,6 +2167,6 @@
        )
        io_group_list(13) % to_be_processed = .FALSE.
-
-       !incoming direct sw flux
+!
+!--    incoming direct sw flux
        io_group_list(14) = init_io_group(                                      &
           in_files = radiation_files,                                          &
@@ -2121,6 +2176,6 @@
        )
        io_group_list(14) % to_be_processed = .FALSE.
-
-       !sw radiation balance
+!
+!--    sw radiation balance
        io_group_list(15) = init_io_group(                                      &
           in_files = radiation_files,                                          &
@@ -2130,6 +2185,6 @@
        )
        io_group_list(15) % to_be_processed = .FALSE.
-
-       !lw radiation balance
+!
+!--    lw radiation balance
        io_group_list(16) = init_io_group(                                      &
           in_files = radiation_files,                                          &
@@ -2143,4 +2198,15 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initializes an IO group with the list of input files, the lists of
+!> input and output variables, and sets the number of output quantities.
+!>
+!> In this context, output quantities refers to the physical quantities required
+!> to interpolate netCDF output variables, not the output variables itsleft. For
+!> instance, the output variables ls_forcing_left/_right/_north/_south all rely
+!> on the output quantity Theta.
+!------------------------------------------------------------------------------!
     FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
                            n_output_quantities) RESULT(group)
@@ -2157,5 +2223,11 @@
        group % n_inputs = SIZE(in_var_list)
        group % kind = TRIM(kind)
-
+!
+!--    For the 'thermodynamics' IO group, one quantity more than input variables
+!--    is needed to compute all output variables of the IO group. Concretely, in
+!--    preprocess() the density is computed from T,P or PP,QV in adddition to
+!--    the variables Theta, p, qv. In read_input_variables(),
+!--    n_output_quantities is used to allocate the correct number of input
+!--    buffers.
        IF ( PRESENT(n_output_quantities) )  THEN
           group % n_output_quantities = n_output_quantities
@@ -2207,5 +2279,5 @@
 ! Description:
 ! ------------
-!> Initializes the the variable list.
+!> Initializes the variable list.
 !------------------------------------------------------------------------------!
     SUBROUTINE setup_variable_tables(ic_mode)
@@ -2787,5 +2859,5 @@
           intermediate_grid = palm_intermediate                                &
        )
-       ! Radiation fluxes and balances
+
        output_var_table(38) = init_nc_var(                                     &
           name              = 'rad_swd_dif_0',                                 &
@@ -3154,5 +3226,6 @@
        output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set
 
-       ! Attributes shared among all variables
+!
+!--    Attributes shared among all variables
        output_var_table(:) % source = nc_source_text
 
@@ -3164,5 +3237,5 @@
 ! Description:
 ! ------------
-!> Initializes and nc_var varible with the given parameters. The 'kind'
+!> Initializes an nc_var varible with the given parameters. The 'kind'
 !> parameter is used to infer the correct netCDF IDs and the level of detail,
 !> 'lod', as defined by the PALM-4U input data standard.
@@ -3199,7 +3272,8 @@
        SELECT CASE( TRIM(out_var_kind) )
 
-       !TODO: Using global module variables 'init_variables_required' and
-       !TODO: 'boundary_variables_required'. Encapsulate in settings type
-       !TODO: and pass into init_nc_var.
+!
+!--    TODO: Using global module variables 'init_variables_required' and
+!--    TODO: 'boundary_variables_required'. Encapsulate in settings type
+!--    TODO: and pass into init_nc_var.
        CASE( 'init soil' )
           var % nt              = 1
@@ -3591,8 +3665,10 @@
        
        CASE( 'velocities' )
-          ! Allocate a compute buffer with the same number of arrays as the input
+!
+!--       Allocate a compute buffer with the same number of arrays as the input
           ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) )
 
-          ! Allocate u and v arrays with scalar dimensions
+!
+!--       Allocate u and v arrays with scalar dimensions
           nx = SIZE(input_buffer(1) % array, 1)
           ny = SIZE(input_buffer(1) % array, 2)
@@ -3603,5 +3679,6 @@
  CALL run_control('time', 'alloc')
 
-          ! interpolate U and V to centres
+!
+!--       interpolate U and V to centres
           CALL centre_velocities( u_face = input_buffer(1) % array,            &
                                   v_face = input_buffer(2) % array,            &
@@ -3613,6 +3690,7 @@
 
           CASE('rotated-pole')
-             ! rotate U and V to PALM-4U orientation and overwrite U and V with
-             ! rotated velocities
+!
+!--          rotate U and V to PALM-4U orientation and overwrite U and V with
+!--          rotated velocities
              DO k = 1, nz
              DO j = 1, ny
@@ -3637,5 +3715,4 @@
           END SELECT
 
-          ! set values 
           input_buffer(1) % array(1,:,:) = 0.0_dp
           input_buffer(2) % array(1,:,:) = 0.0_dp
@@ -3658,5 +3735,6 @@
           nz = SIZE(input_buffer(1) % array, 3)
 
-          ! Compute absolute pressure if presure perturbation has been read in.
+!
+!--       Compute absolute pressure if presure perturbation has been read in.
           IF ( TRIM(group % in_var_list(2) % name) == 'PP' )  THEN
              message = "Absolute pressure, P, not available, " //              &
@@ -3673,4 +3751,7 @@
                                      RD, G, basic_state_pressure)
 
+!
+!--             Overwrite pressure perturbation with absolute pressure. HECTO
+!--             converts pressure perturbation from hPa to Pa.
                 input_buffer (2) % array(i,j,:) =                              &
                    HECTO * input_buffer (2) % array(i,j,:) +                   &
@@ -3687,8 +3768,10 @@
 
           END IF
-          ! pressure
+!
+!--       mark pressure as preprocessed
           input_buffer(2) % is_preprocessed = .TRUE.
 
-          ! Copy temperature to last input buffer array
+!
+!--       Copy temperature to the last input buffer array
           ALLOCATE(                                                            &
               input_buffer( group % n_output_quantities ) % array (nx, ny, nz) &
@@ -3697,5 +3780,6 @@
               input_buffer(1) % array(:,:,:)
 
-          ! Convert absolute to potential temperature
+!
+!--       Convert absolute in place to potential temperature
           CALL potential_temperature(                                          &
              t = input_buffer(1) % array(:,:,:),                               &
@@ -3706,8 +3790,10 @@
           )
 
-          ! potential temperature
+!
+!--       mark potential temperature as preprocessed
           input_buffer(1) % is_preprocessed = .TRUE.
 
-          ! Convert temperature copy to density
+!
+!--       Convert temperature copy to density
           CALL moist_density(                                                  &
              t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), &
@@ -3718,8 +3804,10 @@
           )
 
-          ! qv
+!
+!--       mark qv as preprocessed
           input_buffer(3) % is_preprocessed = .TRUE.
 
-          ! density
+!
+!--       mark density as preprocessed
           input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE.
 
@@ -3780,27 +3868,32 @@
           SELECT CASE(dt)
 
+!
+!--       input has been accumulated over one hour. Leave as is
+!--       input_buffer(1) % array(:,:,:) carrries one-hour integral
           CASE(1)
-             !input has been accumulated over one hour. Leave as is
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral
-
+
+!
+!--       input has been accumulated over two hours. Subtract previous step
+!--       input_buffer(1) % array(:,:,:) carrries one-hour integral
+!--       input_buffer(2) % array(:,:,:) carrries two-hour integral
           CASE(2)
-             !input has been accumulated over two hours. Subtract previous step
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral
-             !input_buffer(2) % array(:,:,:) carrries two-hour integral
              CALL deaverage(                                                   &
                       avg_1 = input_buffer(1) % array(:,:,:), t1 = 1.0_dp,     &
                       avg_2 = input_buffer(2) % array(:,:,:), t2 = 1.0_dp,     &
                       avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp )
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour
-
+!
+!--          input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour
+
+!
+!--       input has been accumulated over three hours. Subtract previous step
+!--       input_buffer(1) % array(:,:,:) carrries three-hour integral
+!--       input_buffer(2) % array(:,:,:) still carrries two-hour integral
           CASE(3)
-             !input has been accumulated over three hours. Subtract previous step
-             !input_buffer(1) % array(:,:,:) carrries three-hour integral
-             !input_buffer(2) % array(:,:,:) still carrries two-hour integral
              CALL deaverage(                                                   &
                      avg_1 = input_buffer(2) % array(:,:,:), t1 = 1.0_dp,      &
                      avg_2 = input_buffer(1) % array(:,:,:), t2 = 1.0_dp,      &
                      avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp )
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA
+!
+!--          input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA
 
           CASE DEFAULT
@@ -3818,28 +3911,34 @@
 
           hour = iter - 1
-          dt = MODULO(hour, 3) + 1 ! averaging period
+!
+!--       averaging period
+          dt = MODULO(hour, 3) + 1
           SELECT CASE(dt)
-
+!
+!--       input has been accumulated over one hour. Leave as is
+!--       input_buffer(1) % array(:,:,:) carrries one-hour integral
           CASE(1)
-             !input has been accumulated over one hour. Leave as is
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral
-
+
+!
+!--       input has been accumulated over two hours. Subtract previous step
+!--       input_buffer(1) % array(:,:,:) carrries one-hour integral
+!--       input_buffer(2) % array(:,:,:) carrries two-hour integral
           CASE(2)
-             !input has been accumulated over two hours. Subtract previous step
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral
-             !input_buffer(2) % array(:,:,:) carrries two-hour integral
              CALL deaverage( input_buffer(1) % array(:,:,:), 1.0_dp,           &
                              input_buffer(2) % array(:,:,:), 2.0_dp,           &
                              input_buffer(1) % array(:,:,:), 1.0_dp)
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour
-
+!
+!--          input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour
+
+!
+!--       input has been accumulated over three hours. Subtract previous step
+!--       input_buffer(1) % array(:,:,:) carrries three-hour integral
+!--       input_buffer(2) % array(:,:,:) still carrries two-hour integral
           CASE(3)
-             !input has been accumulated over three hours. Subtract previous step
-             !input_buffer(1) % array(:,:,:) carrries three-hour integral
-             !input_buffer(2) % array(:,:,:) still carrries two-hour integral
              CALL deaverage( input_buffer(2) % array(:,:,:), 2.0_dp,           &
                              input_buffer(1) % array(:,:,:), 3.0_dp,           &
                              input_buffer(1) % array(:,:,:), 1.0_dp)
-             !input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA
+!
+!--          input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA
 
           CASE DEFAULT
@@ -3860,4 +3959,5 @@
 
 
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -3928,5 +4028,6 @@
                 END DO
 
-                ! Overwrite if at least one non-water neighbour cell is available
+!
+!--             Overwrite if at least one non-water neighbour cell is available
                 IF (n_cells > 0)  THEN
                    array(i,j,:) = column(:) / n_cells
Index: palm/trunk/UTIL/inifor/src/inifor_io.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_io.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor_io.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
+! 
+! 
+! 3537 2018-11-20 10:53:14Z eckhard
 ! New routine get_netcdf_dim_vector()
 ! 
@@ -73,5 +77,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -92,4 +96,12 @@
     IMPLICIT NONE
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
+!> variable 'in_var % name' from the file 'in_file'. The netCDF data array is
+!> stored in the 'buffer' array and metadata added to the respective members of
+!> the given 'in_var'.
+!------------------------------------------------------------------------------!
     INTERFACE get_netcdf_variable
         MODULE PROCEDURE get_netcdf_variable_int
@@ -101,4 +113,10 @@
  CONTAINS
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_variable_int() implements the integer variant for the
+!> get_netcdf_variable interface.
+!------------------------------------------------------------------------------!
     SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
 
@@ -139,4 +157,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_variable_real() implements the real variant for the
+!> get_netcdf_variable interface.
+!------------------------------------------------------------------------------!
     SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
 
@@ -177,9 +201,15 @@
 
 
-    SUBROUTINE get_netcdf_dim_vector(filename, varname, array)
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
+!> netCDF file 'filename'.
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
 
        CHARACTER(LEN=*), INTENT(IN)         ::  filename
-       CHARACTER(LEN=*), INTENT(IN)         ::  varname
-       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  array(:)
+       CHARACTER(LEN=*), INTENT(IN)         ::  coordname
+       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  coords(:)
 
        INTEGER ::  ncid, varid, dimlen
@@ -187,15 +217,15 @@
 
        IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
-            nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR )  THEN
+            nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR )  THEN
 
           CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
           CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
 
-          ALLOCATE(array(dimlen))
-          CALL check(nf90_get_var( ncid, varid, array ))
+          ALLOCATE(coords(dimlen))
+          CALL check(nf90_get_var( ncid, varid, coords))
 
        ELSE
 
-          message = "Failed to read '" // TRIM(varname) // &
+          message = "Failed to read '" // TRIM(coordname) // &
              "' from file '" // TRIM(filename) // "'."
           CALL abort('get_netcdf_dim_vector', message)
@@ -206,4 +236,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
+!> by 'in_var % name' into 'in_var' data structure.
+!------------------------------------------------------------------------------!
     SUBROUTINE get_input_dimensions(in_var, ncid)
 
@@ -232,4 +268,11 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_start_and_count() gets the start position and element counts for
+!> the given netCDF file. This information is used in get_netcdf_variable_int()
+!> and _real() for reading input variables..
+!------------------------------------------------------------------------------!
     SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
 
@@ -251,8 +294,10 @@
        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
+!
+!--       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 reading from second level, e.g. depth = 0.005 instead of 0.0
           start(3) = 2
        END IF
@@ -269,4 +314,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
+!> output.
+!------------------------------------------------------------------------------!
     SUBROUTINE netcdf_define_variable(var, ncid)
 
@@ -286,4 +337,11 @@
     
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
+!> them in the given the 'var' data structure. This information is used later 
+!> for writing output variables in update_output().
+!------------------------------------------------------------------------------!
     SUBROUTINE netcdf_get_dimensions(var, ncid)
 
@@ -293,5 +351,4 @@
         CHARACTER(SNAME)            ::  null
 
-        ! Remember dimension lenghts for NetCDF writing routine
         DO i = 1, var % ndim
            CALL check(nf90_inquire_dimension(ncid, var % dimids(i), &
@@ -306,8 +363,6 @@
 ! Description:
 ! ------------
-!> This routine initializes Inifor. This includes parsing command-line
-!> arguments, setting the names of the input and output file names as well as
-!> the name of the input namelist and, subsequently, reading in and setting grid
-!> parameters for the PALM-4U computational grid.
+!> This routine parses and interpretes the command-line options and stores the
+!> resulting settings in the 'cfg' data structure.
 !------------------------------------------------------------------------------!
     SUBROUTINE parse_command_line_arguments( cfg )
@@ -330,10 +385,4 @@
        IF (arg_count .GT. 0)  THEN
 
-          ! Every option should have an argument.
-          !IF ( MOD(arg_count, 2) .NE. 0 )  THEN
-          !   message = "Syntax error in command line."
-          !   CALL abort('parse_command_line_arguments', message)
-          !END IF
-          
           message = "The -clon and -clat command line options are depricated. " // &
              "Please remove them form your inifor command and specify the " // &
@@ -360,10 +409,8 @@
                 cfg % debug = .TRUE.
 
-             ! Elevation of the PALM-4U domain above sea level
              CASE( '-z0', '-z', '--elevation' )
                 CALL get_option_argument( i, arg )
                 READ(arg, *) cfg % z0
 
-             ! surface pressure, at z0
              CASE( '-p0', '-r', '--surface-pressure' )
                 cfg % p0_is_set = .TRUE.
@@ -371,5 +418,4 @@
                 READ(arg, *) cfg % p0
 
-             ! geostrophic wind in x direction
              CASE( '-ug', '-u', '--geostrophic-u' )
                 cfg % ug_is_set = .TRUE.
@@ -377,5 +423,4 @@
                 READ(arg, *) cfg % ug
 
-             ! geostrophic wind in y direction
              CASE( '-vg', '-v', '--geostrophic-v' )
                 cfg % vg_is_set = .TRUE.
@@ -383,11 +428,6 @@
                 READ(arg, *) cfg % vg
 
-             ! domain centre geographical longitude and latitude
              CASE( '-clon', '-clat' )
                 CALL abort('parse_command_line_arguments', message)          
-                !READ(arg, *) lambda_cg
-                !lambda_cg = lambda_cg * TO_RADIANS
-                !READ(arg, *) phi_cg
-                !phi_cg = phi_cg * TO_RADIANS
 
              CASE( '-path', '-p', '--path' )
@@ -444,10 +484,8 @@
                 cfg % namelist_file = TRIM(arg)
 
-             ! initial condition mode: 'profile' / 'volume'
              CASE( '-mode', '-i', '--init-mode' )
                 CALL get_option_argument( i, arg )
                 cfg % ic_mode = TRIM(arg)
 
-             ! boundary conditions / forcing mode: 'ideal' / 'real'
              CASE( '-f', '--forcing-mode' )
                 CALL get_option_argument( i, arg )
@@ -484,4 +522,10 @@
 
    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Get the argument of the i'th command line option, which is at the location
+!> i+1 of the argument list.
+!------------------------------------------------------------------------------!
    SUBROUTINE get_option_argument(i, arg)
       CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
@@ -494,4 +538,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Checks the INIFOR configuration 'cfg' for plausibility.
+!------------------------------------------------------------------------------!
    SUBROUTINE validate_config(cfg)
       TYPE(inifor_config), INTENT(IN) ::  cfg
@@ -503,5 +552,6 @@
       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.
+!
+!--   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, 'static driver')
@@ -557,4 +607,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check whether the given file is present on the filesystem.
+!------------------------------------------------------------------------------!
    LOGICAL FUNCTION file_present(filename, kind)
       CHARACTER(LEN=PATH), INTENT(IN) ::  filename
@@ -585,9 +640,9 @@
 ! Description:
 ! ------------
-!> This routine initializes the Inifor output file, i.e. the PALM-4U
-!> initializing and forcing data as a NetCDF file.
+!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
+!> file.
 !> 
 !> Besides writing metadata, such as global attributes, coordinates, variables,
-!> in the NetCDF file, various NetCDF IDs are saved for later, when Inifor
+!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
 !> writes the actual data.
 !------------------------------------------------------------------------------!
@@ -611,5 +666,6 @@
        CALL report('setup_netcdf_dimensions', message)
 
-       ! Create the NetCDF file. NF90_CLOBBER selects overwrite mode.
+!
+!--    Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
 #if defined( __netcdf4 )
        CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
@@ -646,7 +702,8 @@
        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)
+!
+!--    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)))
@@ -658,57 +715,79 @@
 !- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
 !------------------------------------------------------------------------------
-       dimids = (/0, 0, 0/) ! reset dimids
-          CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
-          CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
-          CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
-          output_file % dimids_scl = dimids ! save dimids for later
-
-       dimvarids = (/0, 0, 0/) ! reset dimvarids
-          CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
-       output_file % dimvarids_scl = dimvarids ! save dimvarids for later
-
-       ! overwrite third dimid with the one of depth
+!
+!--    reset dimids first
+       dimids = (/0, 0, 0/)
+       CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
+       CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
+       CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
+!
+!--    save dimids for later
+       output_file % dimids_scl = dimids 
+
+!
+!--    reset dimvarids first
+       dimvarids = (/0, 0, 0/)
+       CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
+!
+!--    save dimvarids for later
+       output_file % dimvarids_scl = dimvarids
+
+!
+!--    overwrite third dimid with the one of depth
        CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
-       output_file % dimids_soil = dimids ! save dimids for later
-
-       ! overwrite third dimvarid with the one of depth
+!
+!--    save dimids for later
+       output_file % dimids_soil = dimids
+
+!
+!--    overwrite third dimvarid with the one of depth
        CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
        CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
        CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
        CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
-       output_file % dimvarids_soil = dimvarids ! save dimvarids for later
+!
+!--    save dimvarids for later
+       output_file % dimvarids_soil = dimvarids
 !
 !------------------------------------------------------------------------------
 !- Section 2b: Define dimensions for cell faces/velocities
 !------------------------------------------------------------------------------
-       dimids = (/0, 0, 0/) ! reset dimids
-          CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
-          CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
-          CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
-       output_file % dimids_vel = dimids ! save dimids for later
-
-       dimvarids = (/0, 0, 0/) ! reset dimvarids
-          CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
-       output_file % dimvarids_vel = dimvarids ! save dimvarids for later
+!
+!--    reset dimids first
+       dimids = (/0, 0, 0/)
+       CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
+       CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
+       CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
+!
+!--    save dimids for later
+       output_file % dimids_vel = dimids
+
+!
+!--    reset dimvarids first
+       dimvarids = (/0, 0, 0/)
+       CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
+!
+!--    save dimvarids for later
+       output_file % dimvarids_vel = dimvarids
 
 !
@@ -739,11 +818,14 @@
        CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(3), palm_grid%zw))
 
-       ! TODO Read in soil depths from input file before this.
+!
+!--    TODO Read in soil depths from input file before this.
        CALL check(nf90_put_var(ncid, output_file % dimvarids_soil(3), palm_grid%depths))
 
-       ! Write time vector
+!
+!--    Write time vector
        CALL check(nf90_put_var(ncid, output_file % dimvarid_time, output_file % time))
 
-       ! Close the file
+!
+!--    Close the file
        CALL check(nf90_close(ncid))
 
@@ -751,4 +833,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Defines the netCDF variables to be written to the dynamic driver file
+!------------------------------------------------------------------------------!
     SUBROUTINE setup_netcdf_variables(filename, output_variable_table, debug)
 
@@ -824,6 +911,8 @@
 !- Section 1: Load input buffers for accumulated variables
 !------------------------------------------------------------------------------
+!
+!--    radiation budgets, precipitation
        IF (group % kind == 'running average' .OR.                              &
-           group % kind == 'accumulated')  THEN ! radiation budgets, precipitation
+           group % kind == 'accumulated')  THEN
 
           IF (SIZE(group % in_var_list) .GT. 1 ) THEN
@@ -835,10 +924,14 @@
           END IF
 
-          ! use two buffer arrays
+!
+!--       use two buffer arrays
           nbuffers = 2
           IF ( .NOT. ALLOCATED( buffer ) )  ALLOCATE( buffer(nbuffers) )
 
-          ! chose correct buffer array
-          hour = iter - 1! hour of the day
+!
+!--       hour of the day
+          hour = iter - 1
+!
+!--       chose correct buffer array
           buf_id = select_buffer(hour)
 
@@ -859,20 +952,23 @@
        ELSE
 
-          ! 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)
+!
+!--       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')
           
-          ! Read in all input variables, leave extra buffers-if any-untouched.
+!
+!--       Read in all input variables, leave extra buffers-if any-untouched.
           DO ivar = 1, group % n_inputs
 
              input_var => group % in_var_list(ivar)
 
-             ! Check wheather P or PP is present in input file
+!
+!            Check wheather P or PP is present in input file
              IF (input_var % name == 'P')  THEN
                 input_var % name = TRIM( get_pressure_varname(input_file) )
@@ -891,4 +987,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Select the appropriate buffer ID for accumulated COSMO input variables
+!> depending on the current hour.
+!------------------------------------------------------------------------------!
     INTEGER FUNCTION select_buffer(hour)
        INTEGER, INTENT(IN) ::  hour
@@ -943,4 +1045,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read the given global attribute form the given netCDF file.
+!------------------------------------------------------------------------------!
     FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
 
@@ -948,5 +1055,5 @@
        REAL(dp)                     ::  attribute_value
 
-       INTEGER                      :: ncid
+       INTEGER                      ::  ncid
 
        IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
@@ -966,4 +1073,11 @@
 
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Updates the dynamic driver with the interpolated field of the current
+!> variable at the current time step.
+!------------------------------------------------------------------------------!
     SUBROUTINE update_output(var, array, iter, output_file, cfg)
        TYPE(nc_var), INTENT(IN)  ::  var
@@ -980,5 +1094,6 @@
        )
 
-       ! Skip time dimension for output
+!
+!--    Skip time dimension for output
        ndim = var % ndim
        IF ( var_is_time_dependent )  ndim = var % ndim - 1
@@ -990,5 +1105,6 @@
        CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
 
-       ! Reduce dimension of output array according to variable kind
+!
+!--    Reduce dimension of output array according to variable kind
        SELECT CASE (TRIM(var % kind))
        
@@ -1086,46 +1202,9 @@
 
 
-    SUBROUTINE write_netcdf_variable_2d(var, iter, output_file, buffer)
-       TYPE(nc_var), INTENT(IN)          ::  var
-       INTEGER, INTENT(IN)               ::  iter
-       TYPE(nc_file), INTENT(IN)         ::  output_file
-       REAL(dp), INTENT(IN)              ::  buffer(:,:,:)
-
-       INTEGER ::  ncid, ndim_out, start(4), count(4)
-       LOGICAL ::  last_dimension_is_time
-
-       ndim_out = var % ndim
-
-       last_dimension_is_time = var % dimids( var % ndim ) == output_file % dimid_time
-       IF ( last_dimension_is_time )  THEN
-          ndim_out = ndim_out - 1
-       END IF
-
-       start(:)      = (/1,1,1,iter/)
-       count(1:ndim_out) = var%dimlen(1:ndim_out)
-
-       CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
-
-       IF (TRIM(var % kind) .EQ. 'left/right scalar')  THEN
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer(1,:,:) ) )
-
-       ELSE IF (TRIM(var % kind) .EQ. 'north/south scalar')  THEN
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer(:,1,:) ) )
-
-       ELSE IF (TRIM(var % kind) .EQ. 'top scalar')  THEN
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer(:,:,1) ) )
-       ELSE
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer ) )
-
-       END IF
-       CALL check(nf90_close(ncid))
-
-    END SUBROUTINE write_netcdf_variable_2d
-
-
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Checks the status of a netCDF API call and aborts if an error occured
+!------------------------------------------------------------------------------!
     SUBROUTINE check(status)
 
Index: palm/trunk/UTIL/inifor/src/inifor_transform.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_transform.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor_transform.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+! 
+! 
+! 3537 2018-11-20 10:53:14Z eckhard
 ! bugfix: working precision added
 ! 
@@ -52,5 +56,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -110,8 +114,9 @@
        DO k = nz, LBOUND(out_arr, 3), -1
 
-          ! TODO: Remove IF clause and extrapolate based on a critical vertical
-          ! TODO: index marking the lower bound of COSMO-DE data coverage.
-          ! Check for negative interpolation weights indicating grid points 
-          ! below COSMO-DE domain and extrapolate from the top in such cells.
+!
+!--       TODO: Remove IF clause and extrapolate based on a critical vertical
+!--       TODO: index marking the lower bound of COSMO-DE data coverage.
+!--       Check for negative interpolation weights indicating grid points 
+!--       below COSMO-DE domain and extrapolate from the top in such cells.
           IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz)  THEN
              out_arr(i,j,k) = out_arr(i,j,k+1)
@@ -155,6 +160,7 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar)
-    ! I index 0-based for the indices of the outvar to be consistent with the
-    ! outgrid indices and interpolation weights.
+!
+!--    I index 0-based for the indices of the outvar to be consistent with the
+!--    outgrid indices and interpolation weights.
        TYPE(grid_definition), INTENT(IN)  ::  outgrid
        REAL(dp), INTENT(IN)               ::  invar(0:,0:,0:)
@@ -164,6 +170,7 @@
        INTEGER ::  i, j, k, l
 
-       ! TODO: check if input dimensions are consistent, i.e. ranges are correct
-       IF (UBOUND(outvar, 3) .GT. UBOUND(invar, 3))  THEN
+!
+!--    TODO: check if input dimensions are consistent, i.e. ranges are correct
+       IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) )  THEN
            message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // &
               TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
@@ -190,4 +197,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the horizontal average of the in_arr(:,:,:) and store it in
+!> out_arr(:)
+!------------------------------------------------------------------------------!
     SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
        REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
@@ -200,5 +213,5 @@
        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(' ') // "ii has " // str(SIZE(ii)) // " elements, " //   &
              NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
           CALL abort('average_2d', message)
@@ -211,6 +224,5 @@
              i = ii(l)
              j = jj(l)
-             out_arr(k) = out_arr(k) +&
-                         in_arr(i, j, k)
+             out_arr(k) = out_arr(k) + in_arr(i, j, k)
           END DO
 
@@ -223,4 +235,15 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Three-dimensional interpolation driver. Interpolates from the source_array to
+!> the given palm_grid.
+!>
+!> The routine separates horizontal and vertical
+!> interpolation. In the horizontal interpolation step, the source_array data is
+!> interpolated along COSMO grid levels onto the intermediate grid (vertically
+!> as coarse as COSMO, horizontally as fine as PALM).
+!------------------------------------------------------------------------------!
     SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid)
        TYPE(grid_definition), INTENT(IN) ::  palm_intermediate, palm_grid
@@ -232,15 +255,17 @@
        nx = palm_intermediate % nx
        ny = palm_intermediate % ny
-       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.
+       nlev = palm_intermediate % nz
+
+!
+!--    Interpolate from COSMO to intermediate grid. Allocating with one
+!--    less point in the vertical, since scalars like T have 50 instead of 51
+!--    points in COSMO.
        ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) !
 
        CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
 
-       ! Interpolate from intermediate grid to palm_grid grid, includes
-       ! extrapolation for cells below COSMO-DE domain.
+!
+!--    Interpolate from intermediate grid to palm_grid grid, includes
+!--    extrapolation for cells below COSMO domain.
        CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
 
@@ -250,4 +275,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Average data horizontally from the source_array over the region given by the
+!> averaging grid 'avg_grid' and store the result in 'profile_array'.
+!------------------------------------------------------------------------------!
     SUBROUTINE average_profile(source_array, profile_array, avg_grid)
        TYPE(grid_definition), INTENT(IN)          ::  avg_grid
@@ -265,7 +296,11 @@
           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
+!
+!--       Loop over PALM levels
+          DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1)
+
+!
+!--          Loop over vertical interpolation neighbours
+             DO m = 1, 2
 
                 k_source = avg_grid % kkk(l, k_profile, m)
@@ -274,15 +309,21 @@
                    + 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
+!
+!--          Loop over vertical interpolation neighbours m
+             END DO
+
+!
+!--       Loop over PALM levels k_profile
+          END DO
+
+!
+!--    Loop over horizontal neighbours l
+       END DO
 
        ni_columns = 1.0_dp / avg_grid % n_columns
        profile_array(:) = profile_array(:) * ni_columns
 
-       ! Extrapolate constant to the bottom
+!
+!--    Constant extrapolation to the bottom
        profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
 
@@ -290,4 +331,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Extrapolates density linearly from the level 'k_min' downwards.
+!------------------------------------------------------------------------------!
     SUBROUTINE extrapolate_density(rho, avg_grid)
        REAL(dp), DIMENSION(:), INTENT(INOUT) ::  rho
@@ -308,4 +354,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Driver for extrapolating pressure from PALM level k_min downwards
+!------------------------------------------------------------------------------!
     SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
        REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
@@ -405,5 +456,6 @@
        REAL(dp) :: ri
 
-       ! TODO check dimensions of lat/lon and y/x match
+!
+!--    TODO check dimensions of lat/lon and y/x match
 
        ri = 1.0_dp / r
@@ -438,7 +490,8 @@
        REAL(dp) :: ri
 
-       ! If this elemental function is called with a large array as xy, it is
-       ! computationally more efficient to precompute the inverse radius and
-       ! then muliply.
+!
+!--    If this elemental function is called with a large array as xy, it is
+!--    computationally more efficient to precompute the inverse radius and
+!--    then muliply.
        ri = 1.0_dp / r
 
@@ -448,4 +501,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> For a rotated-pole system with the origin at geographical latitude 'phi_c',
+!> compute the geographical latitude of its rotated north pole.
+!------------------------------------------------------------------------------!
     REAL(dp) FUNCTION phic_to_phin(phi_c)
         REAL(dp), INTENT(IN) ::  phi_c
@@ -456,4 +515,11 @@
 
     
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> For a rotated-pole system with the origin at geographical latitude 'phi_c'
+!> and longitude 'lam_c', compute the geographical longitude of its rotated
+!> north pole.
+!------------------------------------------------------------------------------!
     REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
        REAL(dp), INTENT(IN) ::  phi_c, lam_c
@@ -467,11 +533,21 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set gamma according to whether PALM domain is in the northern or southern
+!> hemisphere of the COSMO rotated-pole system. Gamma assumes either the
+!> value 0 or PI and is needed to work around around a bug in the
+!> rotated-pole coordinate transformations.
+!------------------------------------------------------------------------------!
     REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)
-       REAL(dp), INTENT(IN) ::  phi_cg, phi_ref
-       LOGICAL ::  palm_centre_is_south_of_cosmo_origin
-       
-       palm_centre_is_south_of_cosmo_origin = (phi_cg < phi_ref)
-
-       IF (palm_centre_is_south_of_cosmo_origin)  THEN
+       REAL(dp), INTENT(IN) ::  phi_cg
+       REAL(dp), INTENT(IN) ::  phi_ref
+
+       LOGICAL ::  palm_origin_is_south_of_cosmo_origin
+       
+       palm_origin_is_south_of_cosmo_origin = (phi_cg < phi_ref)
+
+       IF (palm_origin_is_south_of_cosmo_origin)  THEN
            gamma_from_hemisphere = PI
        ELSE
@@ -550,4 +626,9 @@
        
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Rotate the given vector field (x(:), y(:)) by the given 'angle'.
+!------------------------------------------------------------------------------!
     SUBROUTINE rotate_vector_field(x, y, angle)
        REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y  !< x and y coodrinate in arbitrary units
@@ -559,8 +640,9 @@
        sine = SIN(angle * TO_RADIANS)
        cosine = COS(angle * TO_RADIANS)
-       ! RESAHPE() fills columns first, so the rotation matrix becomes
-       !
-       ! rotation = [ cosine   -sine  ]
-       !            [  sine    cosine ]
+!
+!--    RESAHPE() fills columns first, so the rotation matrix becomes
+!--    
+!--    rotation = [ cosine   -sine  ]
+!--               [  sine    cosine ]
        rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
 
@@ -588,5 +670,5 @@
 !>    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)
@@ -664,6 +746,7 @@
        DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
        DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx
-          ! Compute the floating point index corrseponding to PALM-4U grid point
-          ! location along target grid (COSMO-DE) axes.
+!
+!--       Compute the floating point index corrseponding to PALM-4U grid point
+!--       location along target grid (COSMO-DE) axes.
           lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
           latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
@@ -693,4 +776,10 @@
 
     
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computes linear vertical interpolation neighbour indices and weights for each
+!> column of the given palm grid.
+!------------------------------------------------------------------------------!
     SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
                                                             palm_intermediate )
@@ -709,5 +798,6 @@
        nlev = palm_intermediate % nz
 
-       ! in each column of the fine grid, find vertical neighbours of every cell
+!
+!--    in each column of the fine grid, find vertical neighbours of every cell
        DO j = 0, ny
        DO i = 0, nx
@@ -718,12 +808,14 @@
           column_top  = palm_intermediate % h(i,j,nlev)
 
-          ! 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
-          ! the current cell overlaps with the current_height.
+!
+!--       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
+!--       the current cell overlaps with the current_height.
           DO k = 1, nz
 
-             ! Memorize the top and bottom boundaries of the coarse cell and the
-             ! current height within it
+!
+!--          Memorize the top and bottom boundaries of the coarse cell and the
+!--          current height within it
              current_height = palm_grid % z(k) + palm_grid % z0
              h_top    = palm_intermediate % h(i,j,k_intermediate+1)
@@ -738,5 +830,6 @@
              )
 
-             ! set default weights
+!
+!--          set default weights
              palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
 
@@ -755,6 +848,7 @@
 
              ELSE
-                ! cycle through intermediate levels until current
-                ! intermediate-grid cell overlaps with current_height
+!
+!--             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
@@ -777,5 +871,6 @@
                 palm_grid % kk(i,j,k,2) = k_intermediate + 1
 
-                ! copmute vertical weights
+!
+!--             compute vertical weights
                 weight = (h_top - current_height) / (h_top - h_bottom)
                 palm_grid % w_verti(i,j,k,1) = weight
@@ -791,4 +886,16 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computes linear vertical interpolation neighbour indices and weights for each
+!> column of the given averaging grid.
+!>
+!> The difference to the _interp variant of this routine lies in how columns
+!> are adressed. While the _interp variant loops over all PALM grid columns
+!> given by combinations of all index indices (i,j), this routine loops over a
+!> subset of COSMO columns, which are sequantlially stored in the index lists
+!> iii(:) and jjj(:).
+!------------------------------------------------------------------------------!
     SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid )
        TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
@@ -805,5 +912,6 @@
        nlev = SIZE(avg_grid % cosmo_h, 3)
 
-       ! in each column of the fine grid, find vertical neighbours of every cell
+!
+!--    in each column of the fine grid, find vertical neighbours of every cell
        DO l = 1, avg_grid % n_columns
 
@@ -814,18 +922,23 @@
           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.
+!
+!--       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
+!
+!--          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
+!
+!--          COSMO column top is located at 22000m, point_is_above_grid is very
+!--          unlikely.
+             point_is_above_grid = (current_height > column_top)
              point_is_below_grid = (current_height < column_base)
 
@@ -835,5 +948,6 @@
              )
 
-             ! set default weights
+!
+!--          set default weights
              avg_grid % w(l,k_palm,1:2) = 0.0_dp
 
@@ -852,6 +966,7 @@
                 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
+!
+!--             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
@@ -865,6 +980,7 @@
                 END DO
 
-                ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
-                ! k_intermediate = 49 is not the beginning of a valid cell.
+!
+!--             k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
+!--             k_intermediate = 49 is not the beginning of a valid cell.
                 IF (k_intermediate > nlev-1)  THEN
                    message = "Index " // TRIM(str(k_intermediate)) //          &
@@ -876,5 +992,6 @@
                 avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
 
-                ! copmute vertical weights
+!
+!--             compute vertical weights
                 weight = (h_top - current_height) / (h_top - h_bottom)
                 avg_grid % w(l,k_palm,1) = weight
@@ -882,6 +999,11 @@
              END IF
 
-          END DO ! k, PALM levels
-       END DO ! l, averaging columns
+!
+!--       Loop over PALM levels k
+          END DO
+
+!
+!--       Loop over averaging columns l
+       END DO
  
     END SUBROUTINE find_vertical_neighbours_and_weights_average
@@ -931,4 +1053,5 @@
 !                         ii(i,j,1/2)        ii(i,j,3/4)
 !          
+!------------------------------------------------------------------------------!
     SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
        palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
@@ -950,8 +1073,10 @@
        DO i = 0, UBOUND(palm_clon, 1)
       
-          ! weight in lambda direction
+!
+!--       weight in lambda direction
           wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
 
-          ! weight in phi direction
+!
+!--       weight in phi direction
           wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
 
@@ -988,4 +1113,5 @@
 !> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
 !> which means the first centre point has to be omitted and is set to zero.
+!------------------------------------------------------------------------------!
     SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
        REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN)  ::  u_face, v_face
@@ -1004,4 +1130,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the geographical latitude of a point given in rotated-pole cordinates
+!------------------------------------------------------------------------------!
     FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam)
     
@@ -1041,4 +1172,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the geographical latitude of a point given in rotated-pole cordinates
+!------------------------------------------------------------------------------!
     FUNCTION phi2phirot (phi, rla, polphi, pollam)
     
@@ -1072,4 +1208,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the geographical longitude of a point given in rotated-pole cordinates
+!------------------------------------------------------------------------------!
     FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
     
@@ -1123,4 +1264,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the rotated-pole longitude of a point given in geographical cordinates
+!------------------------------------------------------------------------------!
     FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
 
@@ -1162,4 +1308,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Rotate the given velocity vector (u,v) from the geographical to the
+!> rotated-pole system
+!------------------------------------------------------------------------------!
     SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
     
@@ -1187,4 +1339,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Rotate the given velocity vector (urot, vrot) from the rotated-pole to the
+!> geographical system
+!------------------------------------------------------------------------------!
     SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
     
Index: palm/trunk/UTIL/inifor/src/inifor_types.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_types.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor_types.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+! 
+! 
+! 3447 2018-10-29 15:52:54Z eckhard
 ! Renamed source files for compatibilty with PALM build system
 ! 
@@ -51,5 +55,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -66,4 +70,9 @@
  IMPLICIT NONE
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Contaner for the INIFOR command-line configuration
+!------------------------------------------------------------------------------!
  TYPE inifor_config
     CHARACTER(LEN=DATE)  ::  start_date           !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation
@@ -82,27 +91,33 @@
     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
-    CHARACTER(LEN=SNAME) ::  rotation_method
-
-    REAL(dp)             ::  p0
-    REAL(dp)             ::  ug
-    REAL(dp)             ::  vg
-    REAL(dp)             ::  z0                   !< Elevation of the PALM-4U domain above sea level [m]
+    CHARACTER(LEN=SNAME) ::  averaging_mode       !< destinguishes between level-based and heigh-based averaging
+    CHARACTER(LEN=SNAME) ::  bc_mode              !< destinguishes realistic and idealistic forcing
+    CHARACTER(LEN=SNAME) ::  ic_mode              !< destinguishes volume and profile initialization
+    CHARACTER(LEN=SNAME) ::  rotation_method      !< selects method for velocity rotation
+
+    REAL(dp)             ::  p0                   !< manually specified surface pressure [Pa]
+    REAL(dp)             ::  ug                   !< manually spefied geostrophic wind component in x direction [m/s]
+    REAL(dp)             ::  vg                   !< manually spefied geostrophic wind component in y direction [m/s]
+    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  !< 
+    LOGICAL              ::  debug                       !< indicates whether --debug option was given
+    LOGICAL              ::  p0_is_set                   !< indicates whether p0 was set manually
+    LOGICAL              ::  ug_is_set                   !< indicates whether ug was set manually
+    LOGICAL              ::  vg_is_set                   !< indicates whether vg was set manually
+    LOGICAL              ::  flow_prefix_is_set          !<  indicates whether the flow prefix was set manually
+    LOGICAL              ::  input_prefix_is_set         !<  indicates whether the input prefix was set manually
+    LOGICAL              ::  radiation_prefix_is_set     !<  indicates whether the radiation prefix was set manually
+    LOGICAL              ::  soil_prefix_is_set          !<  indicates whether the soil prefix was set manually
+    LOGICAL              ::  soilmoisture_prefix_is_set  !<  indicates whether the soilmoisture prefix was set manually
  END TYPE inifor_config
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Container for grid data, in partucular coordinates, interpolation neighbours
+!> and weights
+!------------------------------------------------------------------------------!
  TYPE grid_definition
     CHARACTER(LEN=SNAME)  ::  name(3)       !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/)
@@ -152,18 +167,28 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Container for name and dimensions of the netCDF output file
+!------------------------------------------------------------------------------!
  TYPE nc_file
-    CHARACTER(LEN=PATH)   ::  name          !< file name
-    INTEGER               ::  dimid_time    !< NetCDF IDs of the time dimension
-    INTEGER               ::  dimids_scl(3) !< NetCDF IDs of the grid dimensions for scalar points x, y, z 
-    INTEGER               ::  dimids_vel(3) !< NetCDF IDs of the grid dimensions for velocity points xu, yu, zu
-    INTEGER               ::  dimids_soil(3)!< NetCDF IDs of the grid dimensions for soil points x, y, depth
-    INTEGER               ::  dimvarid_time !< NetCDF IDs of the time variable
-    INTEGER               ::  dimvarids_scl(3) !< NetCDF IDs of the grid coordinates of scalars x, y, z 
-    INTEGER               ::  dimvarids_vel(3) !< NetCDF IDs of the grid coordinates of velocities xu, yu, zu. Note that velocities are located at mix of both coordinates, e.g. u(xu, y, z).
-    INTEGER               ::  dimvarids_soil(3)!< NetCDF IDs of the grid coordinates for soil points x, y, depth 
-    REAL(dp), POINTER     ::  time(:)       ! vector of output time steps
+    CHARACTER(LEN=PATH)   ::  name              !< file name
+    INTEGER               ::  dimid_time        !< NetCDF IDs of the time dimension
+    INTEGER               ::  dimids_scl(3)     !< NetCDF IDs of the grid dimensions for scalar points x, y, z 
+    INTEGER               ::  dimids_vel(3)     !< NetCDF IDs of the grid dimensions for velocity points xu, yu, zu
+    INTEGER               ::  dimids_soil(3)    !< NetCDF IDs of the grid dimensions for soil points x, y, depth
+    INTEGER               ::  dimvarid_time     !< NetCDF IDs of the time variable
+    INTEGER               ::  dimvarids_scl(3)  !< NetCDF IDs of the grid coordinates of scalars x, y, z 
+    INTEGER               ::  dimvarids_vel(3)  !< NetCDF IDs of the grid coordinates of velocities xu, yu, zu. Note that velocities are located at mix of both coordinates, e.g. u(xu, y, z).
+    INTEGER               ::  dimvarids_soil(3) !< NetCDF IDs of the grid coordinates for soil points x, y, depth 
+    REAL(dp), POINTER     ::  time(:)           !< vector of output time steps
  END TYPE nc_file
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Metadata container for netCDF variables
+!------------------------------------------------------------------------------!
  TYPE nc_var
     INTEGER                               ::  varid     !< NetCDF ID of the variable
@@ -193,5 +218,13 @@
 
 
- 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.
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Input/Output group, groups together nc_var-type output variabels that share 
+!> input variables as well as lists of the netCDF files they are stored in.
+!> For instance, all boundary surfaces and initialization fields of the
+!> potential temperature are base on the input netCDF variables T and P.
+!------------------------------------------------------------------------------!
+ TYPE io_group
     INTEGER                          ::  nt             !< maximum number of output time steps across all output variables
     INTEGER                          ::  nv             !< number of netCDF output variables
@@ -208,7 +241,14 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Container for input data arrays. read_input_variables() allocates a
+!> one-dimensional array of containers, to accomodate all inputs of the given
+!> IO group in one variable.
+!------------------------------------------------------------------------------!
  TYPE container
-   REAL(dp), ALLOCATABLE ::  array(:,:,:)
-   LOGICAL               ::  is_preprocessed = .FALSE.
+   REAL(dp), ALLOCATABLE ::  array(:,:,:)               !< generic data array
+   LOGICAL               ::  is_preprocessed = .FALSE.  !< flag indicating whether input array has been preprocessed
  END TYPE container
 
Index: palm/trunk/UTIL/inifor/src/inifor_util.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_util.f90	(revision 3556)
+++ palm/trunk/UTIL/inifor/src/inifor_util.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+! 
+! 
+! 3447 2018-10-29 15:52:54Z eckhard
 ! Renamed source files for compatibilty with PALM build system
 ! 
@@ -44,5 +48,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -61,4 +65,9 @@
     IMPLICIT NONE
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Fortran implementation of C's struct tm for representing points in time
+!------------------------------------------------------------------------------!
     TYPE, BIND(c) :: tm_struct
        INTEGER(C_INT) :: tm_sec     !< seconds after the minute [0, 61]
@@ -75,4 +84,10 @@
     INTERFACE
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Interface to C's strptime function, which converts a string to a tm
+!> structure.
+!------------------------------------------------------------------------------!
        FUNCTION strptime(string, format, timeinfo) BIND(c, NAME='strptime')
           IMPORT :: C_CHAR, C_SIZE_T, tm_struct
@@ -87,4 +102,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Interface to C's strftime function, which converts the given 'timeinfo'
+!> structure to a string in the given 'format'.
+!------------------------------------------------------------------------------!
        FUNCTION strftime(string, string_len, format, timeinfo) BIND(c, NAME='strftime')
           IMPORT :: C_CHAR, C_SIZE_T, tm_struct
@@ -101,4 +122,12 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Interface to C's mktime function, which converts the given tm structure to
+!> Unix time (number of seconds since 0 UTC, 1st January 1970). INIFOR uses the
+!> side effect that a call to mktime noramlizes the given 'timeinfo' structure,
+!> e.g. increments the date if hours overfow 24.
+!------------------------------------------------------------------------------!
        FUNCTION mktime(timeinfo) BIND(c, NAME='mktime')
           IMPORT :: C_PTR, tm_struct
@@ -115,4 +144,11 @@
  CONTAINS
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Takes a string of the form YYYYMMDDHH, adds the given number of hours to its
+!> tm structure representation, and returns the corresponding string in the same
+!> format
+!------------------------------------------------------------------------------!
     CHARACTER(LEN=DATE) FUNCTION add_hours_to(date_string, hours)
        CHARACTER(LEN=DATE), INTENT(IN)          ::  date_string
@@ -143,4 +179,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Print all members of the given tm structure
+!------------------------------------------------------------------------------!
     SUBROUTINE print_tm(timeinfo)
        TYPE(tm_struct), INTENT(IN) :: timeinfo
@@ -158,4 +199,9 @@
 
     
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize the given tm structure with zero values
+!------------------------------------------------------------------------------!
     SUBROUTINE init_tm(timeinfo)
        TYPE(tm_struct), INTENT(INOUT) :: timeinfo
@@ -177,44 +223,10 @@
 
 
-    SUBROUTINE fake_output_3d(a)
-
-       REAL(dp), INTENT(INOUT)       ::  a(:,:,:)
-       REAL(dp)                      ::  lxi, lyi
-       INTEGER ::  i,j,k
-
-       lyi = 2.0_dp * PI / (SIZE(a, 2) - 1.0_dp)
-       lxi = 2.0_dp * PI / (SIZE(a, 1) - 1.0_dp)
-
-       DO k = 1, SIZE(a, 3)
-       DO j = 1, SIZE(a, 2)
-       DO i = 1, SIZE(a, 1)
-           a(i,j,k) = SIN(lxi * i) * COS(lyi * j) + k
-       END DO
-       END DO
-       END DO
-
-    END SUBROUTINE fake_output_3d
-
-
-    SUBROUTINE fake_output_2d(a, offset)
-
-       REAL(dp), INTENT(INOUT) ::  a(:,:)
-       INTEGER, INTENT(IN)     ::  offset
-       REAL(dp)                ::  lxi, lyi
-       INTEGER                 ::  i,j
-
-       lyi = 2.0_dp*PI / (SIZE(a, 2) - 1.0_dp)
-       lxi = 2.0_dp*PI / (SIZE(a, 1) - 1.0_dp)
-
-       a(:,:) = 1.0_dp
-       DO j = 1, SIZE(a, 2)
-       DO i = 1, SIZE(a, 1)
-          a(i,j) = SIN(lxi * i) * COS(lyi * j) + offset
-       END DO
-       END DO
-
-    END SUBROUTINE fake_output_2d
-
-
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Fill the given array with values equally spaced between and including start
+!> and stop
+!------------------------------------------------------------------------------!
     SUBROUTINE linspace(start, stop, array)
 
@@ -255,4 +267,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> 
+!------------------------------------------------------------------------------!
     SUBROUTINE deaverage(avg_1, t1, avg_2, t2, avg_3, t3)
 
@@ -270,4 +287,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the COSMO-DE/-D2 basic state pressure profile
+!------------------------------------------------------------------------------!
     SUBROUTINE get_basic_state(z, beta, p_sl, t_sl, rd, g, p0)
 
@@ -278,5 +300,5 @@
        REAL(dp), INTENT(IN)  ::  rd     !< ideal gas constant of dry air [J/kg/K]
        REAL(dp), INTENT(IN)  ::  g      !< acceleration of Earth's gravity [m/s^2]
-       REAL(dp), INTENT(OUT) ::  p0(1:) !< COSMO-DE basic state pressure [Pa]
+       REAL(dp), INTENT(OUT) ::  p0(1:) !< COSMO basic state pressure [Pa]
        REAL(dp) ::  root_frac, factor   !< precomputed factors
 
@@ -285,6 +307,6 @@
 
        p0(:) = p_sl * EXP(                                                     &
-                  factor * ( 1.0_dp - SQRT( 1.0_dp - root_frac * z(:) ) )  &
-               )
+                  factor * ( 1.0_dp - SQRT( 1.0_dp - root_frac * z(:) ) )      &
+       )
 
     END SUBROUTINE get_basic_state
@@ -345,4 +367,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Converts the given real value to a string
+!------------------------------------------------------------------------------!
     CHARACTER(LEN=16) FUNCTION real_to_str_f(val)
 
@@ -355,4 +382,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Converts the given integer value to a string
+!------------------------------------------------------------------------------!
     CHARACTER(LEN=10) FUNCTION str(val)
 
@@ -365,17 +397,9 @@
 
 
-    CHARACTER(LEN=30) FUNCTION strs(vals)
-
-        INTEGER, INTENT(IN) ::  vals(:)
-        INTEGER ::  i
-
-        strs = ''
-        DO i = 1, SIZE(vals)
-           strs = TRIM(strs) // TRIM(str(vals(i)))
-        END DO
-
-    END FUNCTION strs
-
-
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> If the given path is not conlcuded by a slash, add one.
+!------------------------------------------------------------------------------!
     SUBROUTINE normalize_path(path)
         
