Index: palm/trunk/UTIL/inifor/src/inifor_grid.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_grid.f90	(revision 3537)
+++ 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
