Changeset 4450 for palm

Timestamp:

Mar 9, 2020 7:12:57 PM (5 years ago)

Author:

suehring

Message:

Bugfix, missing from_file check for a variable from static input file

File:

• palm/trunk/SOURCE/land_surface_model_mod.f90 (modified) (328 diffs)

palm/trunk/SOURCE/land_surface_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Missing from_file check
+!
+! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directive moved
-! 
+!
 ! 4442 2020-03-04 19:21:13Z suehring
-! Change order of dimension in surface arrays %frac, %emissivity and %albedo 
+! Change order of dimension in surface arrays %frac, %emissivity and %albedo
 ! to allow for better vectorization in the radiation interactions.
-! 
+!
 ! 4441 2020-03-04 19:20:35Z suehring
 ! bugfix: missing cpp-directives for serial mode added, misplaced cpp-directives moved
-! 
+!
 ! 4381 2020-01-20 13:51:46Z suehring
-! - Bugfix in nested soil initialization in case no dynamic input file is 
+! - Bugfix in nested soil initialization in case no dynamic input file is
 !   present
-! - In order to do not mess-up the job-protocoll, give error messages 503, 507 
+! - In order to do not mess-up the job-protocoll, give error messages 503, 507
 !   and 508 only once
-! 
+!
 ! 4360 2020-01-07 11:25:50Z suehring
 ! Fix wrong location string in message call
-! 
+!
 ! 4356 2019-12-20 17:09:33Z suehring
 ! Correct single message calls, local checks must be given by the respective
 ! mpi rank.
-! 
+!
 ! 4339 2019-12-13 18:18:30Z suehring
 ! Bugfix, character length too short, caused crash on NEC.
-! 
+!
 ! 4338 2019-12-13 13:23:23Z suehring
-! To avoid divisions by zero, add security factor in calculation of roughness 
+! To avoid divisions by zero, add security factor in calculation of roughness
 ! length over water surfaces.
-! 
+!
 ! 4321 2019-12-04 10:26:38Z pavelkrc
 ! Initialization of relative surface fractions revised
-! 
+!
 ! 4312 2019-11-27 14:06:25Z suehring
 ! Bugfix: partitioning of LE from liquid water reservoir fixed. Bare soils are
 ! now allowed to store liquid water at the surface.
-! 
+!
 ! 4261 2019-10-09 17:58:00Z scharf
 ! bugfix for rev. 4258: deallocate temporary arrays
-! 
+!
 ! 4258 2019-10-07 13:29:08Z suehring
-! - Revise limitation for soil moisture in case it exceeds its saturation 
+! - Revise limitation for soil moisture in case it exceeds its saturation
 !   value (J. Resler)
-! - Revise initialization of soil moisture and temperature in a nested run in 
-!   case dynamic input information is available. This case, the soil within 
+! - Revise initialization of soil moisture and temperature in a nested run in
+!   case dynamic input information is available. This case, the soil within
 !   the child domains can be initialized separately. (J. Resler, M. Suehring)
-! - As part of this revision, migrate the netcdf input of soil temperature / 
-!   moisture to this module, as well as the routine to inter/extrapolate soil 
-!   profiles between different grids. 
-! 
+! - As part of this revision, migrate the netcdf input of soil temperature /
+!   moisture to this module, as well as the routine to inter/extrapolate soil
+!   profiles between different grids.
+!
 ! 4251 2019-10-02 12:07:38Z maronga
 ! Bugfix: albedo_types for vegetation_type look-up table corrected.
-! 
+!
 ! 4201 2019-08-29 15:47:27Z suehring
-! - Limit soil moisture to its saturation moisture and give a respective 
+! - Limit soil moisture to its saturation moisture and give a respective
 !   warning rather than an error.
 ! - Perform checks for soil temperature only when there is no dynamic input
-!   file for the parent or possible child domains. 
-! 
+!   file for the parent or possible child domains.
+!
 ! 4194 2019-08-28 08:09:44Z suehring
-! Apply more strict limitation of z0 over water surfaces in case it exceeds the 
-! surface-layer height, in order to avoid instabilities. 
-! 
+! Apply more strict limitation of z0 over water surfaces in case it exceeds the
+! surface-layer height, in order to avoid instabilities.
+!
 ! 4188 2019-08-26 14:15:47Z suehring
 ! Minor adjustment in error numbers, typos corrected
-! 
+!
 ! 4187 2019-08-26 12:43:15Z suehring
 ! Adjust message call in case of local checks
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4118 2019-07-25 16:11:45Z suehring
 ! Initialization of soil temperature and moisture via dynamic input file only
-! for vegetation and pavement surfaces. 
-! 
+! for vegetation and pavement surfaces.
+!
 ! 4110 2019-07-22 17:05:21Z suehring
-! Relax checks for non-consistent initialization in case static or dynamic 
+! Relax checks for non-consistent initialization in case static or dynamic
 ! input is provided. For example, soil_temperature or deep_soil_temperature
-! is not mandatory any more if dynamic input is available. Also, improper 
-! settings of x_type in namelist are only checked if no static file is 
-! available. 
-! 
+! is not mandatory any more if dynamic input is available. Also, improper
+! settings of x_type in namelist are only checked if no static file is
+! available.
+!
 ! 4109 2019-07-22 17:00:34Z suehring
-! Further revision of last commit in order to avoid any side effects when 
-! albedo type is not set in namelist and default albedo type changes. 
-! 
+! Further revision of last commit in order to avoid any side effects when
+! albedo type is not set in namelist and default albedo type changes.
+!
 ! 4024 2019-06-12 14:06:46Z suehring
 ! Bugfix in albedo initialization, caused crashes in rrtmg calls
-! 
+!
 ! 3987 2019-05-22 09:52:13Z kanani
 ! Introduce alternative switch for debug output during timestepping
-! 
+!
 ! 3964 2019-05-09 09:48:32Z suehring
-! In a nested child domain, distinguish between soil moisture and temperature 
+! In a nested child domain, distinguish between soil moisture and temperature
 ! initialization from parent via dynamic input file. Further, initialize soil
 ! moisture/temperature from dynamic input file only when initialization via
 ! 'inifor' is desired.
-! 
+!
 ! 3943 2019-05-02 09:50:41Z maronga
 ! Removed extra blank character
-! 
+!
 ! 3941 2019-04-30 09:48:33Z suehring
 ! Check that at least one surface type is set at surface element.
-! 
+!
 ! 3933 2019-04-25 12:33:20Z kanani
 ! Remove unused subroutine and allocation of pt_2m, this is done in surface_mod
 ! now (surfaces%pt_2m)
-! 
-! 
-! Changes related to global restructuring of location messages and introduction 
+!
+!
+! Changes related to global restructuring of location messages and introduction
 ! of additional debug messages
-! 
+!
 ! 3881 2019-04-10 09:31:22Z suehring
-! Bugfix in level 3 initialization of pavement albedo type and pavement 
+! Bugfix in level 3 initialization of pavement albedo type and pavement
 ! emissivity
-! 
+!
 ! 3868 2019-04-08 11:52:36Z suehring
-! More strict limitation of roughness length when it is in the order of the 
-! vertical grid spacing 
-! 
+! More strict limitation of roughness length when it is in the order of the
+! vertical grid spacing
+!
 ! 3856 2019-04-03 11:06:59Z suehring
 ! Bugfix in lsm_init in case no surface-fractions are provided
-! 
+!
 ! 3847 2019-04-01 14:51:44Z suehring
-! Adjust message-call for checks that are especially carried out locally. 
-! 
+! Adjust message-call for checks that are especially carried out locally.
+!
 ! 3832 2019-03-28 13:16:58Z raasch
 ! instrumented with openmp directives
-! 
+!
 ! 3786 2019-03-06 16:58:03Z raasch
 ! further unused variables removed
-! 
+!
 ! 3767 2019-02-27 08:18:02Z raasch
 ! unused variable for file index removed from rrd-subroutines parameter list
-! 
+!
 ! 3715 2019-02-04 17:34:55Z suehring
 ! Revise check for saturation moisture
-! 
+!
 ! 3710 2019-01-30 18:11:19Z suehring
 ! Check if soil-, water-, pavement- and vegetation types are set within a valid
 ! range.
-! 
+!
 ! 3692 2019-01-23 14:45:49Z suehring
 ! Revise check for soil moisture higher than its saturation value
-! 
+!
 ! 3685 2019-01-21 01:02:11Z knoop
 ! Some interface calls moved to module_interface + cleanup
-! 
+!
 ! 3677 2019-01-17 09:07:06Z moh.hefny
 ! Removed most_method
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! nopointer option removed
-! 
+!
 ! 1496 2014-12-02 17:25:50Z maronga
 ! Initial revision
-! 
+!
 !
 ! Description:
@@ -188 +191 @@
 !> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
 !> scheme implemented in the ECMWF IFS model, with modifications according to
-!> H-TESSEL. The implementation is based on the formulation implemented in the 
+!> H-TESSEL. The implementation is based on the formulation implemented in the
 !> DALES and UCLA-LES models.
 !>
-!> @todo Extensive verification energy-balance solver for vertical surfaces, 
+!> @todo Extensive verification energy-balance solver for vertical surfaces,
 !>       e.g. parametrization of r_a
-!> @todo Revise single land-surface processes for vertical surfaces, e.g. 
-!>       treatment of humidity, etc. 
-!> @todo Consider partial absorption of the net shortwave radiation by the 
+!> @todo Revise single land-surface processes for vertical surfaces, e.g.
+!>       treatment of humidity, etc.
+!> @todo Consider partial absorption of the net shortwave radiation by the
 !>       skin layer.
 !> @todo Improve surface water parameterization
-!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 
+!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
 !>       nzt_soil=3)).
-!> @todo Implement surface runoff model (required when performing long-term LES 
+!> @todo Implement surface runoff model (required when performing long-term LES
 !>       with considerable precipitation.
 !> @todo Revise calculation of f2 when wilting point is non-constant in the
@@ -207 +210 @@
 !> @note No time step criterion is required as long as the soil layers do not
 !>       become too thin.
-!> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case 
+!> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case
 !>       more levels are used this may cause an potential bug
 !> @todo Routine calc_q_surface required?
@@ -213 +216 @@
 !------------------------------------------------------------------------------!
  MODULE land_surface_model_mod
- 
+
     USE arrays_3d,                                                             &
         ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w, hyrho, exner, d_exner
@@ -235 +238 @@
                surface_pressure, timestep_scheme, tsc,                         &
                time_since_reference_point
-               
+
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point_s
@@ -274 +277 @@
                 vegetation_type_f,                                             &
                 water_pars_f,                                                  &
-                water_type_f             
+                water_type_f
 
     USE kinds
@@ -283 +286 @@
         ONLY:  albedo, albedo_type, emissivity, force_radiation_call,          &
                radiation, radiation_scheme, unscheduled_radiation_calls
-        
+
     USE statistics,                                                            &
         ONLY:  hom, statistic_regions
@@ -303 +306 @@
     REAL(wp), PARAMETER  ::                    &
               b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
-              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil (W/m/K)  
+              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil (W/m/K)
               lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix (W/m/K)
               lambda_h_water     = 0.57_wp,    & ! heat conductivity of water (W/m/K)
@@ -329 +332 @@
 !-- LSM variables
     CHARACTER(10) :: surface_type = 'netcdf'      !< general classification. Allowed are:
-                                                  !< 'vegetation', 'pavement', ('building'), 
+                                                  !< 'vegetation', 'pavement', ('building'),
                                                   !< 'water', and 'netcdf'
 
@@ -339 +342 @@
                     nzs = 8,                  & !< number of soil layers
                     pavement_depth_level = 0, & !< default NAMELIST nzt_pavement
-                    pavement_type = 1,        & !< default NAMELIST pavement_type                 
+                    pavement_type = 1,        & !< default NAMELIST pavement_type
                     soil_type = 3,            & !< default NAMELIST soil_type
                     vegetation_type = 2,      & !< default NAMELIST vegetation_type
                     water_type = 1              !< default NAMELISt water_type
-                    
-   
-        
+
+
+
     LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
                constant_roughness = .FALSE.,     & !< use fixed/dynamic roughness lengths for water surfaces
@@ -390 +393 @@
                 z0_water       = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
                 z0h_water      = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
-                z0q_water      = 9999999.9_wp             !< NAMELIST z0q (lsm_par)  
-                
-                
+                z0q_water      = 9999999.9_wp             !< NAMELIST z0q (lsm_par)
+
+
     REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil_center, & !< 1/dz_soil_center
                                             ddz_soil,        & !< 1/dz_soil
@@ -400 +403 @@
 
 
-                                            
+
     REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp,     & !< (NAMELIST) distribution of root surface area to the individual soil layers
                                    soil_moisture = 0.0_wp,           & !< NAMELIST soil moisture content (m3/m3)
@@ -406 +409 @@
                                    dz_soil  = 9999999.9_wp,          & !< (NAMELIST) soil layer depths (spacing)
                                    zs_layer = 9999999.9_wp             !< soil layer depths (edge)
-                                 
+
     TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
                                      t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
                                      m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
-                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements 
+                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements
 
     TYPE(surf_type_lsm), TARGET  ::  t_soil_h_1,  & !<
@@ -421 +424 @@
                                      t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
                                      m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
-                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements   
+                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements
 
     TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
@@ -429 +432 @@
                                      m_soil_v_2     !<
 
-    TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements 
-                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements 
-                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements 
-                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements 
+    TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
+                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
+                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements
+                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements
 
     TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
@@ -440 +443 @@
 
     TYPE(surf_type_lsm), DIMENSION(:), POINTER  ::    &
-                                      t_surface_v,    & !< surface temperature (K), vertical surface elements 
-                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements 
-                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements 
-                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements 
+                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
+                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
+                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements
+                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements
 
     TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
@@ -456 +459 @@
                                                         m_soil_av    !< Average of m_soil
 
-    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m      !< liquid water reservoir tendency (m), horizontal surface elements 
-    TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements 
-    TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements 
-    TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements 
-
-    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m      !< liquid water reservoir tendency (m), vertical surface elements 
-    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements 
-    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements 
-    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_soil_v_m     !< m_soil storage array, vertical surface elements 
-
-!
-!-- Energy balance variables                
+    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m      !< liquid water reservoir tendency (m), horizontal surface elements
+    TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements
+    TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements
+    TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements
+
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m      !< liquid water reservoir tendency (m), vertical surface elements
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_soil_v_m     !< m_soil storage array, vertical surface elements
+
+!
+!-- Energy balance variables
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
               c_liq_av,         & !< average of c_liq
               c_soil_av,        & !< average of c_soil
               c_veg_av,         & !< average of c_veg
-              lai_av,           & !< average of lai       
+              lai_av,           & !< average of lai
               qsws_liq_av,      & !< average of qsws_liq
               qsws_soil_av,     & !< average of qsws_soil
               qsws_veg_av,      & !< average of qsws_veg
               r_s_av              !< average of r_s
-  
+
 !
 !-- Predefined Land surface classes (vegetation_type)
     CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
-                                   'user defined              ',           & !  0 
-                                   'bare soil                 ',           & !  1                           
+                                   'user defined              ',           & !  0
+                                   'bare soil                 ',           & !  1
                                    'crops, mixed farming      ',           & !  2
                                    'short grass               ',           & !  3
@@ -505 +508 @@
 !-- Soil model classes (soil_type)
     CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
-                                   'user defined',                  & ! 0 
+                                   'user defined',                  & ! 0
                                    'coarse      ',                  & ! 1
                                    'medium      ',                  & ! 2
@@ -517 +520 @@
 !-- Pavement classes
     CHARACTER(29), DIMENSION(0:15), PARAMETER :: pavement_type_name = (/ &
-                                   'user defined                 ', & ! 0 
+                                   'user defined                 ', & ! 0
                                    'asphalt/concrete mix         ', & ! 1
                                    'asphalt (asphalt concrete)   ', & ! 2
@@ -533 +536 @@
                                    'artifical turf (sports)      ', & ! 14
                                    'clay (sports)                '  & ! 15
-                                                                 /)                                                              
-                                                                 
+                                                                 /)
+
 !
 !-- Water classes
     CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
-                                   'user defined',                   & ! 0 
+                                   'user defined',                   & ! 0
                                    'lake        ',                   & ! 1
                                    'river       ',                   & ! 2
@@ -544 +547 @@
                                    'pond        ',                   & ! 4
                                    'fountain    '                    & ! 5
-                                                                  /)                                                                                  
-                    
+                                                                  /)
+
 !
 !-- Land surface parameters according to the respective classes (vegetation_type)
@@ -599 +602 @@
                                                                /), (/ 12, 18 /) )
 
-                                    
+
 !
 !-- Root distribution for default soil layer configuration (sum = 1)
@@ -618 +621 @@
                                  0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 13
                                  0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 14
-                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 15 
-                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16 
+                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 15
+                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 17
                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 18
@@ -628 +631 @@
 
 !
-!-- Soil parameters  alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res 
+!-- Soil parameters  alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res
     REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/     &
                       3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp,& ! 1
@@ -641 +644 @@
 !
 !-- TO BE FILLED
-!-- Pavement parameters      z0,       z0h, albedo_type, emissivity  
+!-- Pavement parameters      z0,       z0h, albedo_type, emissivity
     REAL(wp), DIMENSION(0:3,1:15), PARAMETER :: pavement_pars = RESHAPE( (/ &
                       5.0E-2_wp, 5.0E-4_wp,     18.0_wp,    0.97_wp,  & !  1
                       5.0E-2_wp, 5.0E-4_wp,     19.0_wp,    0.94_wp,  & !  2
-                      1.0E-2_wp, 1.0E-4_wp,     20.0_wp,    0.98_wp,  & !  3                                  
+                      1.0E-2_wp, 1.0E-4_wp,     20.0_wp,    0.98_wp,  & !  3
                       1.0E-2_wp, 1.0E-4_wp,     21.0_wp,    0.93_wp,  & !  4
                       1.0E-2_wp, 1.0E-4_wp,     22.0_wp,    0.97_wp,  & !  5
@@ -651 +654 @@
                       1.0E-2_wp, 1.0E-4_wp,     24.0_wp,    0.97_wp,  & !  7
                       1.0E-2_wp, 1.0E-4_wp,     25.0_wp,    0.94_wp,  & !  8
-                      1.0E-2_wp, 1.0E-4_wp,     26.0_wp,    0.98_wp,  & !  9                                  
+                      1.0E-2_wp, 1.0E-4_wp,     26.0_wp,    0.98_wp,  & !  9
                       1.0E-2_wp, 1.0E-4_wp,     27.0_wp,    0.93_wp,  & ! 10
                       1.0E-2_wp, 1.0E-4_wp,     28.0_wp,    0.97_wp,  & ! 11
@@ -658 +661 @@
                       1.0E-2_wp, 1.0E-4_wp,     31.0_wp,    0.94_wp,  & ! 14
                       1.0E-2_wp, 1.0E-4_wp,     32.0_wp,    0.98_wp   & ! 15
-                      /), (/ 4, 15 /) )                             
+                      /), (/ 4, 15 /) )
 !
 !-- Pavement subsurface parameters part 1: thermal conductivity (W/m/K)
@@ -700 +703 @@
        1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp  & ! 15
                            /), (/ 8, 15 /) )
- 
+
 !
 !-- TO BE FILLED
@@ -710 +713 @@
        283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 4
        283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp  & ! 5
-                                                                     /), (/ 7, 5 /) )                                                                   
-                                                                                                                                      
+                                                                     /), (/ 7, 5 /) )
+
     SAVE
 
@@ -717 +720 @@
     PRIVATE
 
-    
+
 !
 !-- Public functions
     PUBLIC lsm_boundary_condition, lsm_check_data_output,                      &
            lsm_check_data_output_pr,                                           &
-           lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,& 
+           lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,&
            lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
            lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
@@ -745 +748 @@
        MODULE PROCEDURE lsm_check_data_output
     END INTERFACE lsm_check_data_output
-    
+
     INTERFACE lsm_check_data_output_pr
        MODULE PROCEDURE lsm_check_data_output_pr
     END INTERFACE lsm_check_data_output_pr
-    
+
     INTERFACE lsm_check_parameters
        MODULE PROCEDURE lsm_check_parameters
     END INTERFACE lsm_check_parameters
-    
+
     INTERFACE lsm_3d_data_averaging
        MODULE PROCEDURE lsm_3d_data_averaging
@@ -777 +780 @@
        MODULE PROCEDURE lsm_header
     END INTERFACE lsm_header
-    
+
     INTERFACE lsm_init
        MODULE PROCEDURE lsm_init
@@ -785 +788 @@
        MODULE PROCEDURE lsm_init_arrays
     END INTERFACE lsm_init_arrays
-    
+
     INTERFACE lsm_parin
        MODULE PROCEDURE lsm_parin
     END INTERFACE lsm_parin
-    
+
     INTERFACE lsm_soil_model
        MODULE PROCEDURE lsm_soil_model
@@ -812 +815 @@
 ! Description:
 ! ------------
-!> Set internal Neumann boundary condition at outer soil grid points 
-!> for temperature and humidity. 
+!> Set internal Neumann boundary condition at outer soil grid points
+!> for temperature and humidity.
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_boundary_condition
- 
+
     IMPLICIT NONE
 
@@ -879 +882 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
- 
- 
+
+
     USE control_parameters,                                                    &
         ONLY:  data_output, message_string
@@ -886 +889 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  unit  !< 
+    CHARACTER (LEN=*) ::  unit  !<
     CHARACTER (LEN=*) ::  var   !<
 
     INTEGER(iwp) :: i
-    INTEGER(iwp) :: ilen   
+    INTEGER(iwp) :: ilen
     INTEGER(iwp) :: k
 
@@ -902 +905 @@
           ENDIF
           unit = 'm3/m3'
-           
+
        CASE ( 't_soil' )
           IF (  .NOT.  land_surface )  THEN
@@ -909 +912 @@
              CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          unit = 'K'   
-             
+          unit = 'K'
+
        CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'm_liq*',                 &
               'qsws_liq*', 'qsws_soil*', 'qsws_veg*', 'r_s*' )
@@ -969 +972 @@
           ENDIF
 
-          IF ( TRIM( var ) == 'lai*'   )      unit = 'none' 
+          IF ( TRIM( var ) == 'lai*'   )      unit = 'none'
           IF ( TRIM( var ) == 'c_liq*' )      unit = 'none'
           IF ( TRIM( var ) == 'c_soil*')      unit = 'none'
@@ -977 +980 @@
           IF ( TRIM( var ) == 'qsws_soil*' )  unit = 'W/m2'
           IF ( TRIM( var ) == 'qsws_veg*'  )  unit = 'W/m2'
-          IF ( TRIM( var ) == 'r_s*')         unit = 's/m' 
-             
+          IF ( TRIM( var ) == 'r_s*')         unit = 's/m'
+
        CASE DEFAULT
           unit = 'illegal'
@@ -995 +998 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
- 
+
     USE control_parameters,                                                    &
         ONLY:  data_output_pr, message_string
@@ -1006 +1009 @@
 
     IMPLICIT NONE
-   
-    CHARACTER (LEN=*) ::  unit      !< 
-    CHARACTER (LEN=*) ::  variable  !< 
+
+    CHARACTER (LEN=*) ::  unit      !<
+    CHARACTER (LEN=*) ::  variable  !<
     CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
- 
-    INTEGER(iwp) ::  var_count     !< 
+
+    INTEGER(iwp) ::  var_count     !<
 
     SELECT CASE ( TRIM( variable ) )
-       
+
        CASE ( 't_soil', '#t_soil' )
           IF (  .NOT.  land_surface )  THEN
@@ -1059 +1062 @@
 
  END SUBROUTINE lsm_check_data_output_pr
- 
- 
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1070 +1073 @@
     USE control_parameters,                                                    &
         ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string
-                     
-    
+
+
     IMPLICIT NONE
 
@@ -1077 +1080 @@
     INTEGER(iwp) ::  j        !< running index, y-dimension
     INTEGER(iwp) ::  k        !< running index, z-dimension
-    
+
     LOGICAL      ::  dynamic_soil_input_parent !< flag indicating the presence of a dynamic input file for the parent
 
@@ -1086 +1089 @@
          TRIM( surface_type ) /= 'pavement'    .AND.                           &
          TRIM( surface_type ) /= 'water'       .AND.                           &
-         TRIM( surface_type ) /= 'netcdf' )  THEN  
+         TRIM( surface_type ) /= 'netcdf' )  THEN
        message_string = 'unknown surface type: surface_type = "' //            &
                         TRIM( surface_type ) // '"'
@@ -1108 +1111 @@
        CALL message( 'lsm_check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
     ENDIF
-    
+
     IF (  .NOT.  radiation )  THEN
        message_string = 'lsm requires '//                                      &
@@ -1142 +1145 @@
     ENDIF
 !
-!-- Check if vegetation types are set within a valid range.    
+!-- Check if vegetation types are set within a valid range.
     IF ( TRIM( surface_type ) == 'vegetation' )  THEN
        IF ( vegetation_type < LBOUND( vegetation_pars, 2 )  .AND.              &
@@ -1162 +1165 @@
                                         'the valid range at (j,i) = ', j, i
                    CALL message( 'lsm_check_parameters', 'PA0526',             &
-                                  2, 2, myid, 6, 0 ) 
+                                  2, 2, myid, 6, 0 )
                 ENDIF
              ENDDO
@@ -1169 +1172 @@
     ENDIF
 !
-!-- Check if pavement types are set within a valid range.    
+!-- Check if pavement types are set within a valid range.
     IF ( TRIM( surface_type ) == 'pavement' )  THEN
        IF ( pavement_type < LBOUND( pavement_pars, 2 )  .AND.                  &
@@ -1188 +1191 @@
                                         'the valid range at (j,i) = ', j, i
                    CALL message( 'lsm_check_parameters', 'PA0527',             &
-                                  2, 2, myid, 6, 0 ) 
+                                  2, 2, myid, 6, 0 )
                 ENDIF
              ENDDO
@@ -1195 +1198 @@
     ENDIF
 !
-!-- Check if water types are set within a valid range.    
+!-- Check if water types are set within a valid range.
     IF ( TRIM( surface_type ) == 'water' )  THEN
        IF ( water_type < LBOUND( water_pars, 2 )  .AND.                        &
@@ -1214 +1217 @@
                                         'the valid range at (j,i) = ', j, i
                    CALL message( 'lsm_check_parameters', 'PA0528',             &
-                                 2, 2, myid, 6, 0 ) 
+                                 2, 2, myid, 6, 0 )
                 ENDIF
              ENDDO
@@ -1223 +1226 @@
 !-- Check further settings for consistency.
     IF ( TRIM( surface_type ) == 'vegetation' )  THEN
-    
+
        IF ( vegetation_type == 0 )  THEN
           IF ( min_canopy_resistance == 9999999.9_wp )  THEN
@@ -1297 +1300 @@
           ENDIF
        ENDIF
-  
+
     ENDIF
-    
+
     IF ( TRIM( surface_type ) == 'water' )  THEN
 
-       IF ( water_type == 0 )  THEN  
-       
+       IF ( water_type == 0 )  THEN
+
           IF ( z0_water == 9999999.9_wp )  THEN
              message_string = 'water_type = 0 (user_defined)'//                &
@@ -1317 +1320 @@
              CALL message( 'lsm_check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
           ENDIF
-          
+
           IF ( water_temperature == 9999999.9_wp )  THEN
              message_string = 'water_type = 0 (user_defined)'//                &
@@ -1323 +1326 @@
                               '/= 9999999.9'
              CALL message( 'lsm_check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
-          ENDIF       
-          
+          ENDIF
+
        ENDIF
-       
+
     ENDIF
-    
+
     IF ( TRIM( surface_type ) == 'pavement' )  THEN
 
@@ -1337 +1340 @@
        ENDIF
 
-       IF ( pavement_type == 0 )  THEN  
-       
+       IF ( pavement_type == 0 )  THEN
+
           IF ( z0_pavement == 9999999.9_wp )  THEN
              message_string = 'pavement_type = 0 (user_defined)'//             &
@@ -1345 +1348 @@
              CALL message( 'lsm_check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
           ENDIF
-          
+
           IF ( z0h_pavement == 9999999.9_wp )  THEN
              message_string = 'pavement_type = 0 (user_defined)'//             &
@@ -1352 +1355 @@
              CALL message( 'lsm_check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
           ENDIF
-          
+
           IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
              message_string = 'pavement_type = 0 (user_defined)'//             &
@@ -1358 +1361 @@
                               '/= 9999999.9'
              CALL message( 'lsm_check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
-          ENDIF 
-          
+          ENDIF
+
            IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
              message_string = 'pavement_type = 0 (user_defined)'//             &
@@ -1365 +1368 @@
                               '/= 9999999.9'
              CALL message( 'lsm_check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
-          ENDIF 
+          ENDIF
 
           IF ( pavement_depth_level == 0 )  THEN
@@ -1372 +1375 @@
                               '/= 0'
              CALL message( 'lsm_check_parameters', 'PA0474', 1, 2, 0, 6, 0 )
-          ENDIF 
-
-       ENDIF 
-    
+          ENDIF
+
+       ENDIF
+
     ENDIF
 
     IF ( TRIM( surface_type ) == 'netcdf' )  THEN
-!
-!--    MS: Some problme here, after calling message everythings stucks at
-!--        MPI_FINALIZE call. 
-       IF ( ANY( pavement_type_f%var /= pavement_type_f%fill )  .AND.           &
-            ANY( dz_soil /= 9999999.9_wp ) )  THEN
-          message_string = 'pavement-surfaces are not allowed in ' //           &
-                           'combination with a non-default setting of dz_soil'
-          CALL message( 'lsm_check_parameters', 'PA0316', 2, 2, 0, 6, 0 )
+       IF ( pavement_type_f%from_file )  THEN
+          IF ( ANY( pavement_type_f%var /= pavement_type_f%fill )  .AND.       &
+               ANY( dz_soil /= 9999999.9_wp ) )  THEN
+             message_string = 'pavement-surfaces are not allowed in ' //       &
+                              'combination with a non-default setting of dz_soil'
+             CALL message( 'lsm_check_parameters', 'PA0316', 2, 2, 0, 6, 0 )
+          ENDIF
        ENDIF
     ENDIF
-    
+
 !
 !-- Temporary message as long as NetCDF input is not available
@@ -1457 +1459 @@
 
     ENDIF
-
-
-!!! these checks are not needed for water surfaces??
 
 !
@@ -1474 +1473 @@
              nzt_soil = nzt_soil + 1
           ENDIF
-       ENDDO    
+       ENDDO
     ENDIF
     nzs = nzt_soil + 1
 
 !
-!-- Check whether valid soil temperatures are prescribed. Only check this if 
+!-- Check whether valid soil temperatures are prescribed. Only check this if
 !-- no dynamic soil is not initialized with dynamic input.
-!-- In a nested case, check whether there is a dynamic input file for the 
+!-- In a nested case, check whether there is a dynamic input file for the
 !-- child (input_pids_dynamic = .T.) or one for the parent (inquire without
-!-- coupling_char. 
+!-- coupling_char.
     INQUIRE( FILE = TRIM( input_file_dynamic ),                                &
              EXIST = dynamic_soil_input_parent )
@@ -1518 +1517 @@
 !
 !-- Calculate grid spacings. Temperature and moisture are defined at
-!-- the center of the soil layers, whereas gradients/fluxes are 
+!-- the center of the soil layers, whereas gradients/fluxes are
 !-- defined at the edges (_layer)
 !
@@ -1538 +1537 @@
     dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
     zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
- 
+
     DO  k = nzb_soil, nzt_soil-1
        dz_soil_center(k) = zs(k+1) - zs(k)
@@ -1545 +1544 @@
                            '(dz_soil_center(k) <= 0.0)'
           CALL message( 'lsm_check_parameters', 'PA0140', 1, 2, 0, 6, 0 )
-       ENDIF  
+       ENDIF
     ENDDO
-  
+
     dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
-       
+
     ddz_soil_center = 1.0_wp / dz_soil_center
     ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)
@@ -1556 +1555 @@
 
  END SUBROUTINE lsm_check_parameters
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1645 +1644 @@
 
 !
-!-- Index offset of surface element point with respect to adjoining 
-!-- atmospheric grid point 
+!-- Index offset of surface element point with respect to adjoining
+!-- atmospheric grid point
     k_off = surf%koff
     j_off = surf%joff
@@ -1658 +1657 @@
     DO  m = 1, surf%ns
 
-       i   = surf%i(m)            
+       i   = surf%i(m)
        j   = surf%j(m)
        k   = surf%k(m)
@@ -1667 +1666 @@
 !--    parameterization uses a combination of two conductivities: a constant
 !--    conductivity for the skin layer, and a conductivity according to the
-!--    uppermost soil layer. For bare soil and pavements, no skin layer is 
-!--    applied. In these cases, the temperature is assumed to be constant 
-!--    between the surface and the first soil layer. The heat conductivity is 
-!--    then derived from the soil/pavement properties. 
+!--    uppermost soil layer. For bare soil and pavements, no skin layer is
+!--    applied. In these cases, the temperature is assumed to be constant
+!--    between the surface and the first soil layer. The heat conductivity is
+!--    then derived from the soil/pavement properties.
 !--    For water surfaces, the conductivity is already set to 1E10.
 !--    Moreover, the heat capacity is set. For bare soil the heat capacity is
@@ -1684 +1683 @@
           lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
                          lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
-                         
+
           ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) /   &
-                                                     surf%m_sat(nzb_soil,m) ) )                   
-                         
+                                                     surf%m_sat(nzb_soil,m) ) )
+
           lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )   &
                            * ddz_soil(nzb_soil) * 2.0_wp
@@ -1695 +1694 @@
 !--       a heat capacity is that of the soil layer, otherwise it is a
 !--       combination of the conductivities from the skin and the soil layer
-          IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN 
+          IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN
             surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
                               + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
-                              * dz_soil(nzb_soil) * 0.5_wp    
+                              * dz_soil(nzb_soil) * 0.5_wp
             lambda_surface = lambda_soil
 
@@ -1717 +1716 @@
 !--    Set heat capacity of the skin/surface. It is ususally zero when a skin
 !--    layer is used, and non-zero otherwise.
-       c_surface_tmp = surf%c_surface(m) 
+       c_surface_tmp = surf%c_surface(m)
 
 !
@@ -1729 +1728 @@
 !        ELSEIF ( cloud_droplets ) THEN
 !           pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
-!           qv1 = q(k,j,i) 
+!           qv1 = q(k,j,i)
 !        ELSE
 !           pt1 = pt(k,j,i)
@@ -1743 +1742 @@
 !--     heat transfer coefficient for forced convection along vertical walls
 !--     follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
-!--            
+!--
 !--       H = httc (Tsfc - Tair)
 !--       httc = rw * (11.8 + 4.2 * Ueff) - 4.0
-!--            
+!--
 !--             rw: wall patch roughness relative to 1.0 for concrete
 !--             Ueff: effective wind speed
 !--             - 4.0 is a reduction of Rowley et al (1930) formulation based on
 !--             Cole and Sturrock (1977)
-!--           
+!--
 !--             Ucan: Canyon wind speed
 !--             wstar: convective velocity
@@ -1757 +1756 @@
 !--             zH: height of the convective layer
 !--             wstar = (g/Tcan*Qs*zH)**(1./3.)
-                
-!--    Effective velocity components must always 
-!--    be defined at scalar grid point. The wall normal component is 
+
+!--    Effective velocity components must always
+!--    be defined at scalar grid point. The wall normal component is
 !--    obtained by simple linear interpolation. ( An alternative would
 !--    be an logarithmic interpolation. )
@@ -1765 +1764 @@
 !--    1000 is used in the nominator for scaling)
 !--    To do: detailed investigation which approach gives more reliable results!
-!--    Please note, in case of very small friction velocity, e.g. in little 
-!--    holes, the resistance can become negative. For this reason, limit r_a 
-!--    to positive values. 
+!--    Please note, in case of very small friction velocity, e.g. in little
+!--    holes, the resistance can become negative. For this reason, limit r_a
+!--    to positive values.
        IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
           surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
@@ -1778 +1777 @@
                                    ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
                               0.01_wp ) )                                      &
-                           )  - 4.0_wp  ) 
+                           )  - 4.0_wp  )
        ENDIF
 !
-!--    Make sure that the resistance does not drop to zero for neutral 
+!--    Make sure that the resistance does not drop to zero for neutral
 !--    stratification. Also, set a maximum resistance to avoid the breakdown of
 !--    MOST for locations with zero wind speed
        IF ( surf%r_a(m) <   1.0_wp )  surf%r_a(m) =   1.0_wp
-       IF ( surf%r_a(m) > 300.0_wp )  surf%r_a(m) = 300.0_wp       
+       IF ( surf%r_a(m) > 300.0_wp )  surf%r_a(m) = 300.0_wp
 !
 !--    Second step: calculate canopy resistance r_canopy
 !--    f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
- 
-!--    f1: correction for incoming shortwave radiation (stomata close at 
+
+!--    f1: correction for incoming shortwave radiation (stomata close at
 !--    night)
        f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /          &
@@ -1797 +1796 @@
 
 !
-!--    f2: correction for soil moisture availability to plants (the 
+!--    f2: correction for soil moisture availability to plants (the
 !--    integrated soil moisture must thus be considered here)
 !--    f2 = 0 for very dry soils
@@ -1804 +1803 @@
            m_total = m_total + surf%root_fr(ks,m)                              &
                      * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
-       ENDDO 
+       ENDDO
 
 !
@@ -1854 +1853 @@
 
        surf%r_soil(m) = surf%r_soil_min(m) / f2
-       
+
 !
 !--    Calculate the maximum possible liquid water amount on plants and
 !--    bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
-!--    assumed, while paved surfaces might hold up 1 mm of water. The 
-!--    liquid water fraction for paved surfaces is calculated after 
-!--    Masson (2000) (TEB model) and originates from Noilhan & Planton (1989), 
+!--    assumed, while paved surfaces might hold up 1 mm of water. The
+!--    liquid water fraction for paved surfaces is calculated after
+!--    Masson (2000) (TEB model) and originates from Noilhan & Planton (1989),
 !--    while the ECMWF formulation is used for vegetated surfaces and bare soils.
        IF ( surf%pavement_surface(m) )  THEN
@@ -1882 +1881 @@
 
 !
-!--    Calculate coefficients for the total evapotranspiration 
+!--    Calculate coefficients for the total evapotranspiration
 !--    In case of water surface, set vegetation and soil fluxes to zero.
 !--    For pavements, only evaporation of liquid water is possible.
@@ -1914 +1913 @@
 !
 !--    Calculate net radiation radiation without longwave outgoing flux because
-!--    it has a dependency on surface temperature and thus enters the prognostic 
+!--    it has a dependency on surface temperature and thus enters the prognostic
 !--    equations directly
        surf%rad_net_l(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)              &
@@ -1953 +1952 @@
        surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
                           surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
-                                             * dt_3d * tsc(2) ) 
+                                             * dt_3d * tsc(2) )
 
 !
@@ -1985 +1984 @@
 !--    especially when setting skip_time_do_radiation /= 0. The threshold
 !--    value of 0.2 used here is just a first guess. This method should be
-!--    revised in the future as tests have shown that the threshold is 
+!--    revised in the future as tests have shown that the threshold is
 !--    often reached, when no oscillations would occur (causes immense
 !--    computing time for the radiation code).
@@ -2050 +2049 @@
        ENDIF
 !
-!--    Calculate change in liquid water reservoir due to dew fall or 
+!--    Calculate change in liquid water reservoir due to dew fall or
 !--    evaporation of liquid water
        IF ( humidity )  THEN
@@ -2062 +2061 @@
 !--          Add precipitation to liquid water reservoir, if possible.
 !--          Otherwise, add the water to soil. In case of
-!--          pavements, the exceeding water amount is explicitly removed 
+!--          pavements, the exceeding water amount is explicitly removed
 !--          (as fictive runoff by drainage systems)
              IF ( surf%pavement_surface(m) )  THEN
@@ -2070 +2069 @@
                                  * hyrho(k+k_off)                              &
                                  * 0.001_wp * rho_l * l_v
-                ENDIF         
+                ENDIF
              ELSE
                 IF ( surf_m_liq%var_1d(m) < m_liq_max )  THEN
@@ -2080 +2079 @@
                                  surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
                                  * hyrho(k+k_off)                              &
-                                 * 0.001_wp * rho_l * l_v                                 
+                                 * 0.001_wp * rho_l * l_v
                 ELSE
 
@@ -2101 +2100 @@
 !--          Check if reservoir is full (avoid values > m_liq_max)
 !--          In that case, qsws_liq goes to qsws_soil for pervious surfaces. In
-!--          this case qsws_veg is zero anyway (because c_liq = 1),       
+!--          this case qsws_veg is zero anyway (because c_liq = 1),
 !--          so that tend is zero and no further check is needed
              IF ( surf_m_liq%var_1d(m) == m_liq_max )  THEN
@@ -2111 +2110 @@
 
 !
-!--          In case qsws_veg becomes negative (unphysical behavior), 
+!--          In case qsws_veg becomes negative (unphysical behavior),
 !--          let the water enter the liquid water reservoir as dew on the
 !--          plant
@@ -2118 +2117 @@
                 surf%qsws_veg(m) = 0.0_wp
              ENDIF
-          ENDIF                    
- 
+          ENDIF
+
           surf%qsws(m) = surf%qsws(m) / l_v
- 
+
           tend = - surf%qsws_liq(m) * drho_l_lv
           surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *              &
@@ -2176 +2175 @@
 !-- Calculate new roughness lengths (for water surfaces only)
     IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
-    
+
     IF ( debug_output_timestep )  THEN
        WRITE( debug_string, * ) 'lsm_energy_balance', horizontal, l
@@ -2202 +2201 @@
        DO  m = 1, surf%ns
 
-          i   = surf%i(m)            
+          i   = surf%i(m)
           j   = surf%j(m)
           k   = surf%k(m)
 !
 !--       Calculate water vapour pressure at saturation and convert to hPa
-          e_s = 0.01_wp * magnus( MIN(surf_t_surface_p%var_1d(m), 333.15_wp) )                   
+          e_s = 0.01_wp * magnus( MIN(surf_t_surface_p%var_1d(m), 333.15_wp) )
 
 !
@@ -2226 +2225 @@
                                           q(k,j,i)
           ENDIF
-          
+
           surf%q_surface(m) = q(k+k_off,j+j_off,i+i_off)
 !
@@ -2233 +2232 @@
                                   ( 1.0_wp + 0.61_wp * surf%q_surface(m) )
 
-        
-                     
+
+
        ENDDO
        !$OMP END PARALLEL
-  
+
     END SUBROUTINE calc_q_surface
-       
+
  END SUBROUTINE lsm_energy_balance
-    
-   
+
+
 
 !------------------------------------------------------------------------------!
@@ -2260 +2259 @@
        CHARACTER (LEN=86) ::  soil_depth_chr      !< String for soil depth
        CHARACTER (LEN=20) ::  coor_chr            !< Temporary string
-    
+
        INTEGER(iwp) ::  i                         !< Loop index over soil layers
- 
+
        INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
- 
+
        t_soil_chr = ''
        m_soil_chr    = ''
-       soil_depth_chr  = '' 
-       roots_chr        = '' 
+       soil_depth_chr  = ''
+       roots_chr        = ''
        vertical_index_chr   = ''
 
@@ -2313 +2312 @@
 2   FORMAT ('    --> Soil bottom is closed (water content is conserved',       &
             ', default)')
-3   FORMAT ('    --> Soil bottom is open (water content is not conserved)')         
+3   FORMAT ('    --> Soil bottom is open (water content is not conserved)')
 4   FORMAT ('    --> Land surface type  : ',A,/                                &
             '    --> Soil porosity type : ',A)
@@ -2346 +2345 @@
         ONLY:  pmc_is_rootmodel
 #endif
-           
+
        USE pmc_interface,                                                      &
            ONLY:  nested_run
-    
+
        IMPLICIT NONE
 
@@ -2362 +2361 @@
        INTEGER(iwp) ::  m                       !< running index
        INTEGER(iwp) ::  st                      !< soil-type index
-       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4 
+       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4
 #if defined( __parallel )
        INTEGER(iwp) ::  nzs_root                !< number of soil layers in root domain (used in case soil data needs to be
@@ -2374 +2373 @@
 
 #if defined( __parallel )
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  m_soil_root    !< domain-averaged soil moisture profile in root domain 
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  m_soil_root    !< domain-averaged soil moisture profile in root domain
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_soil_root    !< domain-averaged soil temperature profile in root domain
 #endif
@@ -2467 +2466 @@
        ALLOCATE ( surf_lsm_h%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
        ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
-    
+
        surf_lsm_h%lambda_h     = 0.0_wp
 !
@@ -2473 +2472 @@
        IF ( humidity )  THEN
           ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
-          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )  
-
-          surf_lsm_h%lambda_w = 0.0_wp 
+          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
+
+          surf_lsm_h%lambda_w = 0.0_wp
        ENDIF
 !
@@ -2490 +2489 @@
           ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)      )
           ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
-          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )  
-          ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )  
+          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
+          ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
           ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
 
-          surf_lsm_v(l)%lambda_h     = 0.0_wp 
-          
+          surf_lsm_v(l)%lambda_h     = 0.0_wp
+
 !
 !--       If required, allocate humidity-related variables for the soil model
           IF ( humidity )  THEN
              ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
-             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )  
-
-             surf_lsm_v(l)%lambda_w = 0.0_wp 
-          ENDIF     
+             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
+
+             surf_lsm_v(l)%lambda_w = 0.0_wp
+          ENDIF
        ENDDO
 !
-!--    Allocate albedo type and emissivity for vegetation, water and pavement 
+!--    Allocate albedo type and emissivity for vegetation, water and pavement
 !--    fraction.
-!--    Set default values at each surface element. 
+!--    Set default values at each surface element.
        ALLOCATE ( surf_lsm_h%albedo_type(1:surf_lsm_h%ns,0:2) )
        ALLOCATE ( surf_lsm_h%emissivity(1:surf_lsm_h%ns,0:2) )
 !
 !--    Initialize albedo type according to its default type, in order to set values
-!--    independent on default albedo_type in radiation model. 
+!--    independent on default albedo_type in radiation model.
        surf_lsm_h%albedo_type(:,ind_veg_wall)  =                               &
                              INT( vegetation_pars(ind_v_at,vegetation_type) )
@@ -2525 +2524 @@
           ALLOCATE ( surf_lsm_v(l)%emissivity(1:surf_lsm_v(l)%ns,0:2)  )
 !
-!--       Initialize albedo type according to its default type, in order to 
-!--       set values independent on default albedo_type in radiation model. 
+!--       Initialize albedo type according to its default type, in order to
+!--       set values independent on default albedo_type in radiation model.
           surf_lsm_v(l)%albedo_type(:,ind_veg_wall)  =                         &
                              INT( vegetation_pars(ind_v_at,vegetation_type) )
@@ -2536 +2535 @@
        ENDDO
 !
-!--    Allocate arrays for relative surface fraction. 
+!--    Allocate arrays for relative surface fraction.
 !--    0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
        ALLOCATE( surf_lsm_h%frac(1:surf_lsm_h%ns,0:2) )
@@ -2546 +2545 @@
 !
 !--    For vertical walls only - allocate special flag indicating if any building is on
-!--    top of any natural surfaces. Used for initialization only. 
+!--    top of any natural surfaces. Used for initialization only.
        DO  l = 0, 3
           ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) )
@@ -2556 +2555 @@
        ALLOCATE( surf_lsm_h%vegetation_type(1:surf_lsm_h%ns) )
        ALLOCATE( surf_lsm_h%water_type(1:surf_lsm_h%ns)      )
-       
+
        surf_lsm_h%pavement_type   = 0
        surf_lsm_h%vegetation_type = 0
        surf_lsm_h%water_type      = 0
-       
+
        ALLOCATE( surf_lsm_h%pavement_type_name(1:surf_lsm_h%ns)   )
        ALLOCATE( surf_lsm_h%vegetation_type_name(1:surf_lsm_h%ns) )
        ALLOCATE( surf_lsm_h%water_type_name(1:surf_lsm_h%ns)      )
-       
+
        surf_lsm_h%pavement_type_name   = 'none'
        surf_lsm_h%vegetation_type_name = 'none'
        surf_lsm_h%water_type_name      = 'none'
-       
+
        DO  l = 0, 3
           ALLOCATE( surf_lsm_v(l)%pavement_type(1:surf_lsm_v(l)%ns)   )
           ALLOCATE( surf_lsm_v(l)%vegetation_type(1:surf_lsm_v(l)%ns) )
           ALLOCATE( surf_lsm_v(l)%water_type(1:surf_lsm_v(l)%ns)      )
-          
+
           surf_lsm_v(l)%pavement_type   = 0
           surf_lsm_v(l)%vegetation_type = 0
           surf_lsm_v(l)%water_type      = 0
-       
+
           ALLOCATE( surf_lsm_v(l)%pavement_type_name(1:surf_lsm_v(l)%ns)   )
           ALLOCATE( surf_lsm_v(l)%vegetation_type_name(1:surf_lsm_v(l)%ns) )
           ALLOCATE( surf_lsm_v(l)%water_type_name(1:surf_lsm_v(l)%ns)      )
-       
+
           surf_lsm_v(l)%pavement_type_name   = 'none'
           surf_lsm_v(l)%vegetation_type_name = 'none'
-          surf_lsm_v(l)%water_type_name      = 'none'       
+          surf_lsm_v(l)%water_type_name      = 'none'
        ENDDO
-       
-!
-!--    Set flag parameter for the prescribed surface type depending on user 
+
+!
+!--    Set flag parameter for the prescribed surface type depending on user
 !--    input. Set surface fraction to 1 for the respective type.
        SELECT CASE ( TRIM( surface_type ) )
-          
+
           CASE ( 'vegetation' )
-         
+
              surf_lsm_h%vegetation_surface = .TRUE.
              surf_lsm_h%frac(:,ind_veg_wall) = 1.0_wp
@@ -2600 +2599 @@
                 surf_lsm_v(l)%frac(:,ind_veg_wall) = 1.0_wp
              ENDDO
-   
+
           CASE ( 'water' )
-             
+
              surf_lsm_h%water_surface = .TRUE.
              surf_lsm_h%frac(:,ind_wat_win) = 1.0_wp
 !
 !--          Note, vertical water surface does not really make sense.
-             DO  l = 0, 3  
+             DO  l = 0, 3
                 surf_lsm_v(l)%water_surface   = .TRUE.
                 surf_lsm_v(l)%frac(:,ind_wat_win) = 1.0_wp
@@ -2613 +2612 @@
 
           CASE ( 'pavement' )
-             
+
              surf_lsm_h%pavement_surface = .TRUE.
                 surf_lsm_h%frac(:,ind_pav_green) = 1.0_wp
@@ -2644 +2643 @@
                                   2, 2, myid, 6, 0 )
                 ENDIF
-                
+
              ENDDO
 !
 !--          For vertical surfaces some special checks and treatment are
-!--          required for correct initialization. 
+!--          required for correct initialization.
              DO  l = 0, 3
                 DO  m = 1, surf_lsm_v(l)%ns
 !
-!--                Only for vertical surfaces. Check if at the grid point where 
-!--                the wall is defined (i+ioff, j+joff) is any building. 
+!--                Only for vertical surfaces. Check if at the grid point where
+!--                the wall is defined (i+ioff, j+joff) is any building.
 !--                This case, no natural surfaces properties will be defined at
 !--                at this grid point, leading to problems in the initialization.
 !--                To overcome this, define a special flag which
-!--                indicates that a building is defined at the wall grid point  
-!--                and take the surface properties from the adjoining grid 
+!--                indicates that a building is defined at the wall grid point
+!--                and take the surface properties from the adjoining grid
 !--                point, i.e. without offset values.
 !--                Further, there can occur a special case where elevation
 !--                changes are larger than building heights. This case, (j,i)
-!--                and (j+joff,i+ioff) grid points may be both covered by 
-!--                buildings, but vertical, but vertically natural walls may 
-!--                be located between the buildings. This case, it is not 
-!--                guaranteed that information about natural surface types is 
-!--                given, neither at (j,i) nor at (j+joff,i+ioff), again leading 
-!--                to non-initialized surface properties. 
+!--                and (j+joff,i+ioff) grid points may be both covered by
+!--                buildings, but vertical, but vertically natural walls may
+!--                be located between the buildings. This case, it is not
+!--                guaranteed that information about natural surface types is
+!--                given, neither at (j,i) nor at (j+joff,i+ioff), again leading
+!--                to non-initialized surface properties.
                    surf_lsm_v(l)%building_covered(m) = .FALSE.
 !
-!--                Wall grid point is building-covered. This case, set 
-!--                flag indicating that surface properties are initialized 
-!--                from neighboring reference grid point, which is not 
-!--                building_covered. 
+!--                Wall grid point is building-covered. This case, set
+!--                flag indicating that surface properties are initialized
+!--                from neighboring reference grid point, which is not
+!--                building_covered.
                    IF ( building_type_f%from_file )  THEN
                       i = surf_lsm_v(l)%i(m)
@@ -2685 +2684 @@
 !--                   point are both building-covered. This case, surface
 !--                   properties are not necessarily defined (not covered by
-!--                   checks for static input file) at this surface. Hence, 
-!--                   initialize surface properties by simply setting 
-!--                   vegetation_type_f to bare-soil bulk parametrization. 
-!--                   soil_type_f as well as surface_fractions_f will be set 
-!--                   also.                      
+!--                   checks for static input file) at this surface. Hence,
+!--                   initialize surface properties by simply setting
+!--                   vegetation_type_f to bare-soil bulk parametrization.
+!--                   soil_type_f as well as surface_fractions_f will be set
+!--                   also.
                       IF ( building_type_f%var(j+surf_lsm_v(l)%joff,           &
                                                i+surf_lsm_v(l)%ioff) /=        &
@@ -2696 +2695 @@
                       THEN
                          vegetation_type_f%var(j,i)                 = 1 ! bare soil
-                         soil_type_f%var_2d(j,i)                    = 1 
-!                         
-!--                      If surface_fraction is provided in static input, 
+                         soil_type_f%var_2d(j,i)                    = 1
+!
+!--                      If surface_fraction is provided in static input,
 !--                      set fraction for vegetation to one at building-covered
 !--                      surfaces.
@@ -2707 +2706 @@
                          ENDIF
                       ENDIF
-                      
+
                    ENDIF
 !
@@ -2738 +2737 @@
        END SELECT
 !
-!--    In case of netcdf input file, further initialize surface fractions. 
+!--    In case of netcdf input file, further initialize surface fractions.
 !--    At the moment only 1 surface is given at a location, so that the fraction
-!--    is either 0 or 1. This will be revised later. If surface fraction 
+!--    is either 0 or 1. This will be revised later. If surface fraction
 !--    is not given in static input file, relative fractions will be derived
 !--    from given surface type. In this case, only 1 type is given at a certain
-!--    location (already checked).  
+!--    location (already checked).
        IF ( input_pids_static  .AND.  surface_fraction_f%from_file )  THEN
           DO  m = 1, surf_lsm_h%ns
@@ -2773 +2772 @@
                                  j, i, ') are all zero.'
                 CALL message( 'land_surface_model_mod', 'PA0688',              &
-                               2, 2, myid, 6, 0 ) 
+                               2, 2, myid, 6, 0 )
              ENDIF
 !
 !--          In case the sum of all surfaces is not 1, which may happen
 !--          due to rounding errors or type conversions, normalize the
-!--          fractions to one. Note, at the moment no tile approach is 
+!--          fractions to one. Note, at the moment no tile approach is
 !--          implemented, so that relative fractions are either 1 or zero.
              IF ( SUM ( surf_lsm_h%frac(m,:) ) > 1.0_wp  .OR.                  &
@@ -2791 +2790 @@
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-!
-!--             0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction        
+                                                surf_lsm_v(l)%building_covered(m) )
+!
+!--             0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction
                 IF ( surface_fraction_f%frac(ind_veg_wall,j,i) /=              &
                      surface_fraction_f%fill )  THEN
@@ -2819 +2818 @@
                                  j, i, ') are all zero.'
                    CALL message( 'land_surface_model_mod', 'PA0688',           &
-                                  2, 2, myid, 6, 0 ) 
+                                  2, 2, myid, 6, 0 )
                 ENDIF
 !
 !--             In case the sum of all surfaces is not 1, which may happen
 !--             due to rounding errors or type conversions, normalize the
-!--             fractions to one. Note, at the moment no tile approach is 
+!--             fractions to one. Note, at the moment no tile approach is
 !--             implemented, so that relative fractions are either 1 or zero.
                 IF ( SUM ( surf_lsm_v(l)%frac(m,:) ) > 1.0_wp  .OR.            &
@@ -2840 +2839 @@
              j = surf_lsm_h%j(m)
 
-             IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )       &        
+             IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )       &
                 surf_lsm_h%frac(m,ind_veg_wall)  = 1.0_wp
-             IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )       &        
-                surf_lsm_h%frac(m,ind_pav_green) = 1.0_wp 
-             IF ( water_type_f%var(j,i)      /= water_type_f%fill      )       &        
-                surf_lsm_h%frac(m,ind_wat_win)   = 1.0_wp       
+             IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )       &
+                surf_lsm_h%frac(m,ind_pav_green) = 1.0_wp
+             IF ( water_type_f%var(j,i)      /= water_type_f%fill      )       &
+                surf_lsm_h%frac(m,ind_wat_win)   = 1.0_wp
           ENDDO
           DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-     
-                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &        
+                                                surf_lsm_v(l)%building_covered(m) )
+
+                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &
                    surf_lsm_v(l)%frac(m,ind_veg_wall)  = 1.0_wp
-                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )    &        
-                   surf_lsm_v(l)%frac(m,ind_pav_green) = 1.0_wp 
-                IF ( water_type_f%var(j,i)      /= water_type_f%fill      )    &        
-                   surf_lsm_v(l)%frac(m,ind_wat_win)   = 1.0_wp     
+                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )    &
+                   surf_lsm_v(l)%frac(m,ind_pav_green) = 1.0_wp
+                IF ( water_type_f%var(j,i)      /= water_type_f%fill      )    &
+                   surf_lsm_v(l)%frac(m,ind_wat_win)   = 1.0_wp
              ENDDO
           ENDDO
@@ -2865 +2864 @@
 !
 !--    Level 1, initialization of soil parameters.
-!--    It is possible to overwrite each parameter by setting the respecticy 
-!--    NAMELIST variable to a value /= 9999999.9. 
-       IF ( soil_type /= 0 )  THEN  
- 
+!--    It is possible to overwrite each parameter by setting the respecticy
+!--    NAMELIST variable to a value /= 9999999.9.
+       IF ( soil_type /= 0 )  THEN
+
           IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
              alpha_vangenuchten = soil_pars(0,soil_type)
@@ -2878 +2877 @@
 
           IF ( n_vangenuchten == 9999999.9_wp )  THEN
-             n_vangenuchten = soil_pars(2,soil_type)            
+             n_vangenuchten = soil_pars(2,soil_type)
           ENDIF
 
           IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
-             hydraulic_conductivity = soil_pars(3,soil_type)            
+             hydraulic_conductivity = soil_pars(3,soil_type)
           ENDIF
 
           IF ( saturation_moisture == 9999999.9_wp )  THEN
-             saturation_moisture = soil_pars(4,soil_type)           
+             saturation_moisture = soil_pars(4,soil_type)
           ENDIF
 
           IF ( field_capacity == 9999999.9_wp )  THEN
-             field_capacity = soil_pars(5,soil_type)           
+             field_capacity = soil_pars(5,soil_type)
           ENDIF
 
           IF ( wilting_point == 9999999.9_wp )  THEN
-             wilting_point = soil_pars(6,soil_type)            
+             wilting_point = soil_pars(6,soil_type)
           ENDIF
 
           IF ( residual_moisture == 9999999.9_wp )  THEN
-             residual_moisture = soil_pars(7,soil_type)       
+             residual_moisture = soil_pars(7,soil_type)
           ENDIF
 
@@ -2907 +2906 @@
        surf_lsm_h%alpha_vg      = alpha_vangenuchten
        surf_lsm_h%l_vg          = l_vangenuchten
-       surf_lsm_h%n_vg          = n_vangenuchten 
+       surf_lsm_h%n_vg          = n_vangenuchten
        surf_lsm_h%gamma_w_sat   = hydraulic_conductivity
        surf_lsm_h%m_sat         = saturation_moisture
@@ -2919 +2918 @@
           surf_lsm_v(l)%alpha_vg      = alpha_vangenuchten
           surf_lsm_v(l)%l_vg          = l_vangenuchten
-          surf_lsm_v(l)%n_vg          = n_vangenuchten 
+          surf_lsm_v(l)%n_vg          = n_vangenuchten
           surf_lsm_v(l)%gamma_w_sat   = hydraulic_conductivity
           surf_lsm_v(l)%m_sat         = saturation_moisture
@@ -2929 +2928 @@
 !
 !--    Level 2, initialization of soil parameters via soil_type read from file.
-!--    Soil parameters are initialized for each (y,x)-grid point 
+!--    Soil parameters are initialized for each (y,x)-grid point
 !--    individually using default paramter settings according to the given
 !--    soil type.
        IF ( soil_type_f%from_file )  THEN
 !
-!--       Level of detail = 1, i.e. a homogeneous soil distribution along the 
+!--       Level of detail = 1, i.e. a homogeneous soil distribution along the
 !--       vertical dimension is assumed.
           IF ( soil_type_f%lod == 1 )  THEN
@@ -2942 +2941 @@
                 i = surf_lsm_h%i(m)
                 j = surf_lsm_h%j(m)
-             
+
                 st = soil_type_f%var_2d(j,i)
                 IF ( st /= soil_type_f%fill )  THEN
@@ -2956 +2955 @@
              ENDDO
 !
-!--          Vertical surfaces ( assumes the soil type given at respective (x,y) 
+!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
              DO  l = 0, 3
                 DO  m = 1, surf_lsm_v(l)%ns
                    i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
                    j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
 
                    st = soil_type_f%var_2d(j,i)
@@ -2979 +2978 @@
 !
 !--       Level of detail = 2, i.e. soil type and thus the soil parameters
-!--       can be heterogeneous along the vertical dimension. 
+!--       can be heterogeneous along the vertical dimension.
           ELSE
 !
@@ -2986 +2985 @@
                 i = surf_lsm_h%i(m)
                 j = surf_lsm_h%j(m)
-             
+
                 DO  k = nzb_soil, nzt_soil
                    st = soil_type_f%var_3d(k,j,i)
@@ -3002 +3001 @@
              ENDDO
 !
-!--          Vertical surfaces ( assumes the soil type given at respective (x,y) 
+!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
              DO  l = 0, 3
                 DO  m = 1, surf_lsm_v(l)%ns
                    i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
                    j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
 
                    DO  k = nzb_soil, nzt_soil
@@ -3028 +3027 @@
        ENDIF
 !
-!--    Level 3, initialization of single soil parameters at single z,x,y 
+!--    Level 3, initialization of single soil parameters at single z,x,y
 !--    position via soil_pars read from file.
        IF ( soil_pars_f%from_file )  THEN
 !
-!--       Level of detail = 1, i.e. a homogeneous vertical distribution of soil 
+!--       Level of detail = 1, i.e. a homogeneous vertical distribution of soil
 !--       parameters is assumed.
 !--       Horizontal surfaces
           IF ( soil_pars_f%lod == 1 )  THEN
 !
-!--          Horizontal surfaces 
+!--          Horizontal surfaces
              DO  m = 1, surf_lsm_h%ns
                 i = surf_lsm_h%i(m)
@@ -3061 +3060 @@
              ENDDO
 !
-!--          Vertical surfaces 
+!--          Vertical surfaces
              DO  l = 0, 3
                 DO  m = 1, surf_lsm_v(l)%ns
                    i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
                    j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
 
                    IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )           &
@@ -3089 +3088 @@
              ENDDO
 !
-!--       Level of detail = 2, i.e. soil parameters can be set at each soil 
-!--       layer individually. 
+!--       Level of detail = 2, i.e. soil parameters can be set at each soil
+!--       layer individually.
           ELSE
 !
@@ -3119 +3118 @@
              ENDDO
 !
-!--          Vertical surfaces 
+!--          Vertical surfaces
              DO  l = 0, 3
                 DO  m = 1, surf_lsm_v(l)%ns
                    i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
                    j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
 
                    DO  k = nzb_soil, nzt_soil
@@ -3153 +3152 @@
 
 !
-!--    Level 1, initialization of vegetation parameters. A horizontally 
-!--    homogeneous distribution is assumed here. 
+!--    Level 1, initialization of vegetation parameters. A horizontally
+!--    homogeneous distribution is assumed here.
        IF ( vegetation_type /= 0 )  THEN
 
@@ -3162 +3161 @@
 
           IF ( leaf_area_index == 9999999.9_wp )  THEN
-             leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)         
+             leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)
           ENDIF
 
           IF ( vegetation_coverage == 9999999.9_wp )  THEN
-             vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)     
+             vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)
           ENDIF
 
           IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
-              canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)     
+              canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)
           ENDIF
 
           IF ( z0_vegetation == 9999999.9_wp )  THEN
-             z0_vegetation  = vegetation_pars(ind_v_z0,vegetation_type) 
+             z0_vegetation  = vegetation_pars(ind_v_z0,vegetation_type)
           ENDIF
 
           IF ( z0h_vegetation == 9999999.9_wp )  THEN
              z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
-          ENDIF 
-          
+          ENDIF
+
           IF ( z0q_vegetation == 9999999.9_wp )  THEN
              z0q_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
           ENDIF
-          
+
           IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-             lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type) 
+             lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type)
           ENDIF
 
           IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-             lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)           
+             lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)
           ENDIF
 
           IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
-             f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)        
+             f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)
           ENDIF
 
           IF ( c_surface == 9999999.9_wp )  THEN
-             c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)        
+             c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)
           ENDIF
 
           IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
-             albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))       
-          ENDIF
-   
+             albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))
+          ENDIF
+
           IF ( emissivity == 9999999.9_wp )  THEN
-             emissivity = vegetation_pars(ind_v_emis,vegetation_type)     
+             emissivity = vegetation_pars(ind_v_emis,vegetation_type)
           ENDIF
 
@@ -3227 +3226 @@
              surf_lsm_h%albedo_type(m,ind_veg_wall) = albedo_type
              surf_lsm_h%emissivity(m,ind_veg_wall)  = emissivity
-             
+
              surf_lsm_h%vegetation_type(m)      = vegetation_type
              surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
@@ -3235 +3234 @@
              surf_lsm_h%g_d(m)   = 0.0_wp
           ENDIF
-  
+
        ENDDO
 !
@@ -3256 +3255 @@
                 surf_lsm_v(l)%albedo_type(m,ind_veg_wall) = albedo_type
                 surf_lsm_v(l)%emissivity(m,ind_veg_wall)  = emissivity
-                
+
                 surf_lsm_v(l)%vegetation_type(m)      = vegetation_type
                 surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
@@ -3269 +3268 @@
 !
 !--    Level 2, initialization of vegation parameters via vegetation_type read
-!--    from file. Vegetation parameters are initialized for each (y,x)-grid point 
+!--    from file. Vegetation parameters are initialized for each (y,x)-grid point
 !--    individually using default paramter settings according to the given
 !--    vegetation type.
@@ -3278 +3277 @@
              i = surf_lsm_h%i(m)
              j = surf_lsm_h%j(m)
-             
+
              st = vegetation_type_f%var(j,i)
              IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
@@ -3294 +3293 @@
                 surf_lsm_h%albedo_type(m,ind_veg_wall) = INT( vegetation_pars(ind_v_at,st) )
                 surf_lsm_h%emissivity(m,ind_veg_wall)  = vegetation_pars(ind_v_emis,st)
-                
+
                 surf_lsm_h%vegetation_type(m)      = st
                 surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(st)
@@ -3304 +3303 @@
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-             
+                                                surf_lsm_v(l)%building_covered(m) )
+
                 st = vegetation_type_f%var(j,i)
                 IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
@@ -3323 +3322 @@
                    surf_lsm_v(l)%albedo_type(m,ind_veg_wall) = INT( vegetation_pars(ind_v_at,st) )
                    surf_lsm_v(l)%emissivity(m,ind_veg_wall)  = vegetation_pars(ind_v_emis,st)
-                   
+
                    surf_lsm_v(l)%vegetation_type(m)      = st
                    surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(st)
@@ -3331 +3330 @@
        ENDIF
 !
-!--    Level 3, initialization of vegation parameters at single (x,y) 
+!--    Level 3, initialization of vegation parameters at single (x,y)
 !--    position via vegetation_pars read from file.
        IF ( vegetation_pars_f%from_file )  THEN
 !
-!--       Horizontal surfaces 
+!--       Horizontal surfaces
           DO  m = 1, surf_lsm_h%ns
 
@@ -3341 +3340 @@
              j = surf_lsm_h%j(m)
 !
-!--          If surface element is not a vegetation surface and any value in 
+!--          If surface element is not a vegetation surface and any value in
 !--          vegetation_pars is given, neglect this information and give an
-!--          informative message that this value will not be used.   
+!--          informative message that this value will not be used.
              IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.                &
                    ANY( vegetation_pars_f%pars_xy(:,j,i) /=                    &
@@ -3351 +3350 @@
                                  j, i, ') is not a vegetation surface, ',      &
                                  'so that information given in ',              &
-                                 'vegetation_pars at this point is neglected.' 
+                                 'vegetation_pars at this point is neglected.'
                 CALL message( 'land_surface_model_mod', 'PA0436', 0, 0, myid, 6, 0 )
              ELSE
@@ -3409 +3408 @@
           ENDDO
 !
-!--       Vertical surfaces 
+!--       Vertical surfaces
           DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-!
-!--             If surface element is not a vegetation surface and any value in 
+                                                surf_lsm_v(l)%building_covered(m) )
+!
+!--             If surface element is not a vegetation surface and any value in
 !--             vegetation_pars is given, neglect this information and give an
-!--             informative message that this value will not be used.   
+!--             informative message that this value will not be used.
                 IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.          &
                       ANY( vegetation_pars_f%pars_xy(:,j,i) /=                 &
@@ -3427 +3426 @@
                                  j, i, ') is not a vegetation surface, ',      &
                                  'so that information given in ',              &
-                                 'vegetation_pars at this point is neglected.' 
+                                 'vegetation_pars at this point is neglected.'
                    CALL message( 'land_surface_model_mod', 'PA0436', 0, 0, myid, 6, 0 )
                 ELSE
@@ -3486 +3485 @@
              ENDDO
           ENDDO
-       ENDIF 
-
-!
-!--    Level 1, initialization of water parameters. A horizontally 
-!--    homogeneous distribution is assumed here. 
+       ENDIF
+
+!
+!--    Level 1, initialization of water parameters. A horizontally
+!--    homogeneous distribution is assumed here.
        IF ( water_type /= 0 )  THEN
 
           IF ( water_temperature == 9999999.9_wp )  THEN
-             water_temperature = water_pars(ind_w_temp,water_type)       
+             water_temperature = water_pars(ind_w_temp,water_type)
           ENDIF
 
           IF ( z0_water == 9999999.9_wp )  THEN
-             z0_water = water_pars(ind_w_z0,water_type)       
-          ENDIF       
+             z0_water = water_pars(ind_w_z0,water_type)
+          ENDIF
 
           IF ( z0h_water == 9999999.9_wp )  THEN
-             z0h_water = water_pars(ind_w_z0h,water_type)       
-          ENDIF  
-          
+             z0h_water = water_pars(ind_w_z0h,water_type)
+          ENDIF
+
           IF ( z0q_water == 9999999.9_wp )  THEN
-             z0q_water = water_pars(ind_w_z0h,water_type)       
-          ENDIF 
+             z0q_water = water_pars(ind_w_z0h,water_type)
+          ENDIF
 
           IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
-             albedo_type = INT(water_pars(ind_w_at,water_type))        
-          ENDIF
-   
+             albedo_type = INT(water_pars(ind_w_at,water_type))
+          ENDIF
+
           IF ( emissivity == 9999999.9_wp )  THEN
-             emissivity = water_pars(ind_w_emis,water_type)        
-          ENDIF
-
-       ENDIF 
+             emissivity = water_pars(ind_w_emis,water_type)
+          ENDIF
+
+       ENDIF
 !
 !--    Map values onto horizontal elemements
@@ -3528 +3527 @@
              surf_lsm_h%z0q(m)              = z0q_water
              surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
-             surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
+             surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp
              surf_lsm_h%c_surface(m)        = 0.0_wp
              surf_lsm_h%albedo_type(m,ind_wat_win) = albedo_type
              surf_lsm_h%emissivity(m,ind_wat_win)  = emissivity
-             
+
              surf_lsm_h%water_type(m)      = water_type
              surf_lsm_h%water_type_name(m) = water_type_name(water_type)
@@ -3549 +3548 @@
                 surf_lsm_v(l)%z0q(m)              = z0q_water
                 surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
-                surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp               
+                surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp
                 surf_lsm_v(l)%c_surface(m)        = 0.0_wp
                 surf_lsm_v(l)%albedo_type(m,ind_wat_win) = albedo_type
                 surf_lsm_v(l)%emissivity(m,ind_wat_win)  = emissivity
-                
+
                 surf_lsm_v(l)%water_type(m)      = water_type
                 surf_lsm_v(l)%water_type_name(m) = water_type_name(water_type)
-             ENDIF 
+             ENDIF
           ENDDO
        ENDDO
@@ -3562 +3561 @@
 !
 !--    Level 2, initialization of water parameters via water_type read
-!--    from file. Water surfaces are initialized for each (y,x)-grid point 
+!--    from file. Water surfaces are initialized for each (y,x)-grid point
 !--    individually using default paramter settings according to the given
 !--    water type.
@@ -3573 +3572 @@
              i = surf_lsm_h%i(m)
              j = surf_lsm_h%j(m)
-             
+
              st = water_type_f%var(j,i)
              IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
@@ -3582 +3581 @@
                 surf_lsm_h%z0q(m)    = water_pars(ind_w_z0h,st)
                 surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st)
-                surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)              
+                surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)
                 surf_lsm_h%c_surface(m)        = 0.0_wp
                 surf_lsm_h%albedo_type(m,ind_wat_win) = INT( water_pars(ind_w_at,st) )
                 surf_lsm_h%emissivity(m,ind_wat_win)  = water_pars(ind_w_emis,st)
-                
+
                 surf_lsm_h%water_type(m)      = st
                 surf_lsm_h%water_type_name(m) = water_type_name(st)
@@ -3596 +3595 @@
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-             
+                                                surf_lsm_v(l)%building_covered(m) )
+
                 st = water_type_f%var(j,i)
                 IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
@@ -3610 +3609 @@
                                                    water_pars(ind_w_lambda_s,st)
                    surf_lsm_v(l)%lambda_surface_u(m) =                         &
-                                                   water_pars(ind_w_lambda_u,st)           
+                                                   water_pars(ind_w_lambda_u,st)
                    surf_lsm_v(l)%c_surface(m)     = 0.0_wp
                    surf_lsm_v(l)%albedo_type(m,ind_wat_win) =                  &
@@ -3616 +3615 @@
                    surf_lsm_v(l)%emissivity(m,ind_wat_win)  =                  &
                                                   water_pars(ind_w_emis,st)
-                                                  
+
                    surf_lsm_v(l)%water_type(m)      = st
                    surf_lsm_v(l)%water_type_name(m) = water_type_name(st)
@@ -3622 +3621 @@
              ENDDO
           ENDDO
-       ENDIF      
-
-!
-!--    Level 3, initialization of water parameters at single (x,y) 
+       ENDIF
+
+!
+!--    Level 3, initialization of water parameters at single (x,y)
 !--    position via water_pars read from file.
        IF ( water_pars_f%from_file )  THEN
 !
-!--       Horizontal surfaces 
+!--       Horizontal surfaces
           DO  m = 1, surf_lsm_h%ns
              i = surf_lsm_h%i(m)
              j = surf_lsm_h%j(m)
 !
-!--          If surface element is not a water surface and any value in 
+!--          If surface element is not a water surface and any value in
 !--          water_pars is given, neglect this information and give an
-!--          informative message that this value will not be used.   
+!--          informative message that this value will not be used.
              IF ( .NOT. surf_lsm_h%water_surface(m)  .AND.                     &
                    ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) )  THEN
@@ -3643 +3642 @@
                               j, i, ') is not a water surface, ',              &
                               'so that information given in ',                 &
-                              'water_pars at this point is neglected.' 
+                              'water_pars at this point is neglected.'
                 CALL message( 'land_surface_model_mod', 'PA0645', 0, 0, myid, 6, 0 )
              ELSE
@@ -3667 +3666 @@
                       water_pars_f%fill )                                      &
                    surf_lsm_h%lambda_surface_u(m) =                            &
-                                        water_pars_f%pars_xy(ind_w_lambda_u,j,i)     
-       
+                                        water_pars_f%pars_xy(ind_w_lambda_u,j,i)
+
                 IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                     &
                      water_pars_f%fill )                                       &
@@ -3677 +3676 @@
                      water_pars_f%fill )                                       &
                    surf_lsm_h%emissivity(m,ind_wat_win) =                      &
-                                          water_pars_f%pars_xy(ind_w_emis,j,i)  
+                                          water_pars_f%pars_xy(ind_w_emis,j,i)
              ENDIF
           ENDDO
 !
-!--       Vertical surfaces 
+!--       Vertical surfaces
           DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-!
-!--             If surface element is not a water surface and any value in 
+                                                surf_lsm_v(l)%building_covered(m) )
+!
+!--             If surface element is not a water surface and any value in
 !--             water_pars is given, neglect this information and give an
-!--             informative message that this value will not be used.   
+!--             informative message that this value will not be used.
                 IF ( .NOT. surf_lsm_v(l)%water_surface(m)  .AND.               &
                       ANY( water_pars_f%pars_xy(:,j,i) /=                      &
@@ -3699 +3698 @@
                               j, i, ') is not a water surface, ',              &
                               'so that information given in ',                 &
-                              'water_pars at this point is neglected.' 
+                              'water_pars at this point is neglected.'
                    CALL message( 'land_surface_model_mod', 'PA0645',           &
                                   0, 0, myid, 6, 0 )
@@ -3727 +3726 @@
                         water_pars_f%fill )                                    &
                       surf_lsm_v(l)%lambda_surface_u(m) =                      &
-                                      water_pars_f%pars_xy(ind_w_lambda_u,j,i)    
- 
+                                      water_pars_f%pars_xy(ind_w_lambda_u,j,i)
+
                    IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                  &
                         water_pars_f%fill )                                    &
@@ -3737 +3736 @@
                         water_pars_f%fill )                                    &
                       surf_lsm_v(l)%emissivity(m,ind_wat_win)  =               &
-                                      water_pars_f%pars_xy(ind_w_emis,j,i)  
+                                      water_pars_f%pars_xy(ind_w_emis,j,i)
                 ENDIF
              ENDDO
@@ -3745 +3744 @@
 !
 !--    Initialize pavement-type surfaces, level 1
-       IF ( pavement_type /= 0 )  THEN  
-
-!
-!--       When a pavement_type is used, overwrite a possible setting of 
+       IF ( pavement_type /= 0 )  THEN
+
+!
+!--       When a pavement_type is used, overwrite a possible setting of
 !--       the pavement depth as it is already defined by the pavement type
           pavement_depth_level = 0
 
           IF ( z0_pavement == 9999999.9_wp )  THEN
-             z0_pavement  = pavement_pars(ind_p_z0,pavement_type) 
+             z0_pavement  = pavement_pars(ind_p_z0,pavement_type)
           ENDIF
 
@@ -3759 +3758 @@
              z0h_pavement = pavement_pars(ind_p_z0h,pavement_type)
           ENDIF
-          
+
           IF ( z0q_pavement == 9999999.9_wp )  THEN
              z0q_pavement = pavement_pars(ind_p_z0h,pavement_type)
@@ -3770 +3769 @@
           IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
              pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
-          ENDIF   
-    
+          ENDIF
+
           IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
-             albedo_type = INT(pavement_pars(ind_p_at,pavement_type))        
-          ENDIF
-   
+             albedo_type = INT(pavement_pars(ind_p_at,pavement_type))
+          ENDIF
+
           IF ( emissivity == 9999999.9_wp )  THEN
-             emissivity = pavement_pars(ind_p_emis,pavement_type)        
+             emissivity = pavement_pars(ind_p_emis,pavement_type)
           ENDIF
 
@@ -3784 +3783 @@
 !--       lookup table.
           IF ( pavement_depth_level == 0 )  THEN
-             DO  k = nzb_soil, nzt_soil  
+             DO  k = nzb_soil, nzt_soil
                 IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
                 .OR. pavement_subsurface_pars_2(k,pavement_type) == 9999999.9_wp)&
@@ -3798 +3797 @@
        ENDIF
 !
-!--    Level 1 initialization of pavement type surfaces. Horizontally 
+!--    Level 1 initialization of pavement type surfaces. Horizontally
 !--    homogeneous characteristics are assumed
        surf_lsm_h%nzt_pavement = pavement_depth_level
@@ -3809 +3808 @@
              surf_lsm_h%lambda_surface_s(m)    = pavement_heat_conduct         &
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp    
+                                                  * 2.0_wp
              surf_lsm_h%lambda_surface_u(m)    = pavement_heat_conduct         &
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp           
+                                                  * 2.0_wp
              surf_lsm_h%c_surface(m)           = pavement_heat_capacity        &
                                                         * dz_soil(nzb_soil)    &
-                                                        * 0.25_wp                                   
+                                                        * 0.25_wp
 
              surf_lsm_h%albedo_type(m,ind_pav_green) = albedo_type
-             surf_lsm_h%emissivity(m,ind_pav_green)  = emissivity      
-             
+             surf_lsm_h%emissivity(m,ind_pav_green)  = emissivity
+
              surf_lsm_h%pavement_type(m)      = pavement_type
              surf_lsm_h%pavement_type_name(m) = pavement_type_name(pavement_type)
-      
+
              IF ( pavement_type /= 0 )  THEN
                 DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
                    surf_lsm_h%lambda_h_def(k,m)    =                           &
-                                     pavement_subsurface_pars_1(k,pavement_type)                       
+                                     pavement_subsurface_pars_1(k,pavement_type)
                    surf_lsm_h%rho_c_total_def(k,m) =                           &
-                                     pavement_subsurface_pars_2(k,pavement_type) 
+                                     pavement_subsurface_pars_2(k,pavement_type)
                 ENDDO
              ELSE
                 surf_lsm_h%lambda_h_def(:,m)     = pavement_heat_conduct
                 surf_lsm_h%rho_c_total_def(:,m)  = pavement_heat_capacity
-             ENDIF        
-          ENDIF
-       ENDDO                               
+             ENDIF
+          ENDIF
+       ENDDO
 
        DO  l = 0, 3
@@ -3847 +3846 @@
                 surf_lsm_v(l)%lambda_surface_s(m)    = pavement_heat_conduct   &
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp    
+                                                  * 2.0_wp
                 surf_lsm_v(l)%lambda_surface_u(m)    = pavement_heat_conduct   &
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp           
+                                                  * 2.0_wp
                 surf_lsm_v(l)%c_surface(m)           = pavement_heat_capacity  &
                                                         * dz_soil(nzb_soil)    &
-                                                        * 0.25_wp                                      
+                                                        * 0.25_wp
 
                 surf_lsm_v(l)%albedo_type(m,ind_pav_green) = albedo_type
                 surf_lsm_v(l)%emissivity(m,ind_pav_green)  = emissivity
-                
+
                 surf_lsm_v(l)%pavement_type(m)      = pavement_type
                 surf_lsm_v(l)%pavement_type_name(m) = pavement_type_name(pavement_type)
@@ -3864 +3863 @@
                    DO  k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m)
                       surf_lsm_v(l)%lambda_h_def(k,m)    =                     &
-                                     pavement_subsurface_pars_1(k,pavement_type)                       
+                                     pavement_subsurface_pars_1(k,pavement_type)
                       surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
-                                     pavement_subsurface_pars_2(k,pavement_type) 
+                                     pavement_subsurface_pars_2(k,pavement_type)
                    ENDDO
                 ELSE
                    surf_lsm_v(l)%lambda_h_def(:,m)     = pavement_heat_conduct
                    surf_lsm_v(l)%rho_c_total_def(:,m)  = pavement_heat_capacity
-                ENDIF     
+                ENDIF
              ENDIF
           ENDDO
@@ -3877 +3876 @@
 !
 !--    Level 2 initialization of pavement type surfaces via pavement_type read
-!--    from file. Pavement surfaces are initialized for each (y,x)-grid point 
+!--    from file. Pavement surfaces are initialized for each (y,x)-grid point
 !--    individually.
        IF ( pavement_type_f%from_file )  THEN
@@ -3885 +3884 @@
              i = surf_lsm_h%i(m)
              j = surf_lsm_h%j(m)
-             
+
              st = pavement_type_f%var(j,i)
              IF ( st /= pavement_type_f%fill  .AND.  st /= 0 )  THEN
 !
 !--             Determine deepmost index of pavement layer
-                DO  k = nzb_soil, nzt_soil  
+                DO  k = nzb_soil, nzt_soil
                    IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp       &
                    .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)      &
@@ -3906 +3905 @@
                                               pavement_subsurface_pars_1(0,st) &
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp    
+                                                  * 2.0_wp
                 surf_lsm_h%lambda_surface_u(m)  =                              &
                                               pavement_subsurface_pars_1(0,st) &
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp        
+                                                  * 2.0_wp
                 surf_lsm_h%c_surface(m)         =                              &
                                                pavement_subsurface_pars_2(0,st)&
                                                         * dz_soil(nzb_soil)    &
-                                                        * 0.25_wp                               
+                                                        * 0.25_wp
                 surf_lsm_h%albedo_type(m,ind_pav_green) = INT( pavement_pars(ind_p_at,st) )
-                surf_lsm_h%emissivity(m,ind_pav_green)  = pavement_pars(ind_p_emis,st)  
-                
+                surf_lsm_h%emissivity(m,ind_pav_green)  = pavement_pars(ind_p_emis,st)
+
                 surf_lsm_h%pavement_type(m)      = st
                 surf_lsm_h%pavement_type_name(m) = pavement_type_name(st)
@@ -3923 +3922 @@
                 DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
                    surf_lsm_h%lambda_h_def(k,m)    =                           &
-                                     pavement_subsurface_pars_1(k,pavement_type)                       
+                                     pavement_subsurface_pars_1(k,pavement_type)
                    surf_lsm_h%rho_c_total_def(k,m) =                           &
-                                     pavement_subsurface_pars_2(k,pavement_type) 
-                ENDDO   
+                                     pavement_subsurface_pars_2(k,pavement_type)
+                ENDDO
              ENDIF
           ENDDO
@@ -3934 +3933 @@
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-             
+                                                surf_lsm_v(l)%building_covered(m) )
+
                 st = pavement_type_f%var(j,i)
                 IF ( st /= pavement_type_f%fill  .AND.  st /= 0 )  THEN
 !
 !--                Determine deepmost index of pavement layer
-                   DO  k = nzb_soil, nzt_soil  
+                   DO  k = nzb_soil, nzt_soil
                       IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp    &
                       .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)   &
@@ -3957 +3956 @@
                    surf_lsm_v(l)%lambda_surface_s(m)  =                        &
                                               pavement_subsurface_pars_1(0,st) &
-                                                  * ddz_soil(nzb_soil)         & 
-                                                  * 2.0_wp    
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp
                    surf_lsm_v(l)%lambda_surface_u(m)  =                        &
                                               pavement_subsurface_pars_1(0,st) &
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp     
+                                                  * 2.0_wp
 
                    surf_lsm_v(l)%c_surface(m)    =                             &
                                            pavement_subsurface_pars_2(0,st)    &
                                                         * dz_soil(nzb_soil)    &
-                                                        * 0.25_wp                                    
+                                                        * 0.25_wp
                    surf_lsm_v(l)%albedo_type(m,ind_pav_green) =                &
                                               INT( pavement_pars(ind_p_at,st) )
                    surf_lsm_v(l)%emissivity(m,ind_pav_green)  =                &
-                                              pavement_pars(ind_p_emis,st)  
-                                              
+                                              pavement_pars(ind_p_emis,st)
+
                    surf_lsm_v(l)%pavement_type(m)      = st
                    surf_lsm_v(l)%pavement_type_name(m) = pavement_type_name(st)
-                                              
+
                    DO  k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m)
                       surf_lsm_v(l)%lambda_h_def(k,m)    =                     &
-                                    pavement_subsurface_pars_1(k,pavement_type)                       
+                                    pavement_subsurface_pars_1(k,pavement_type)
                       surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
-                                    pavement_subsurface_pars_2(k,pavement_type) 
-                   ENDDO   
+                                    pavement_subsurface_pars_2(k,pavement_type)
+                   ENDDO
                 ENDIF
              ENDDO
           ENDDO
-       ENDIF 
-!
-!--    Level 3, initialization of pavement parameters at single (x,y) 
+       ENDIF
+!
+!--    Level 3, initialization of pavement parameters at single (x,y)
 !--    position via pavement_pars read from file.
        IF ( pavement_pars_f%from_file )  THEN
 !
-!--       Horizontal surfaces 
+!--       Horizontal surfaces
           DO  m = 1, surf_lsm_h%ns
              i = surf_lsm_h%i(m)
              j = surf_lsm_h%j(m)
 !
-!--          If surface element is not a pavement surface and any value in 
+!--          If surface element is not a pavement surface and any value in
 !--          pavement_pars is given, neglect this information and give an
-!--          informative message that this value will not be used.   
+!--          informative message that this value will not be used.
              IF ( .NOT. surf_lsm_h%pavement_surface(m)  .AND.                  &
                    ANY( pavement_pars_f%pars_xy(:,j,i) /=                      &
@@ -4006 +4005 @@
                               j, i, ') is not a pavement surface, ',           &
                               'so that information given in ',                 &
-                              'pavement_pars at this point is neglected.' 
+                              'pavement_pars at this point is neglected.'
                 CALL message( 'land_surface_model_mod', 'PA0647', 0, 0, myid, 6, 0 )
              ELSE
@@ -4026 +4025 @@
                       pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp    
+                                                  * 2.0_wp
                 ENDIF
                 IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) /= &
@@ -4033 +4032 @@
                       pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)   &
                                                   * dz_soil(nzb_soil)          &
-                                                  * 0.25_wp                                    
+                                                  * 0.25_wp
                 ENDIF
                 IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=                  &
@@ -4043 +4042 @@
                    surf_lsm_h%emissivity(m,ind_pav_green)  =                   &
                                    pavement_pars_f%pars_xy(ind_p_emis,j,i)
-             ENDIF 
+             ENDIF
 
           ENDDO
 !
-!--       Vertical surfaces 
+!--       Vertical surfaces
           DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
                 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
-                                                surf_lsm_v(l)%building_covered(m) ) 
-!
-!--             If surface element is not a pavement surface and any value in 
+                                                surf_lsm_v(l)%building_covered(m) )
+!
+!--             If surface element is not a pavement surface and any value in
 !--             pavement_pars is given, neglect this information and give an
-!--             informative message that this value will not be used.   
+!--             informative message that this value will not be used.
                 IF ( .NOT. surf_lsm_v(l)%pavement_surface(m)  .AND.            &
                       ANY( pavement_pars_f%pars_xy(:,j,i) /=                   &
@@ -4065 +4064 @@
                                  j, i, ') is not a pavement surface, ',        &
                                  'so that information given in ',              &
-                                 'pavement_pars at this point is neglected.' 
+                                 'pavement_pars at this point is neglected.'
                    CALL message( 'land_surface_model_mod', 'PA0647', 0, 0, myid, 6, 0 )
                 ELSE
@@ -4086 +4085 @@
                       pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
                                                   * ddz_soil(nzb_soil)         &
-                                                  * 2.0_wp    
+                                                  * 2.0_wp
                    ENDIF
                    IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) &
@@ -4093 +4092 @@
                          pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)&
                                                   * dz_soil(nzb_soil)          &
-                                                  * 0.25_wp                                  
+                                                  * 0.25_wp
                    ENDIF
                    IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=               &
@@ -4103 +4102 @@
                         pavement_pars_f%fill )                                 &
                       surf_lsm_v(l)%emissivity(m,ind_pav_green)  =             &
-                                   pavement_pars_f%pars_xy(ind_p_emis,j,i)  
-                ENDIF 
+                                   pavement_pars_f%pars_xy(ind_p_emis,j,i)
+                ENDIF
              ENDDO
           ENDDO
        ENDIF
 !
-!--    Moreover, for grid points which are flagged with pavement-type 0 or whre 
-!--    pavement_subsurface_pars_f is provided, soil heat conductivity and 
-!--    capacity are initialized with parameters given in        
-!--    pavement_subsurface_pars read from file. 
+!--    Moreover, for grid points which are flagged with pavement-type 0 or whre
+!--    pavement_subsurface_pars_f is provided, soil heat conductivity and
+!--    capacity are initialized with parameters given in
+!--    pavement_subsurface_pars read from file.
        IF ( pavement_subsurface_pars_f%from_file )  THEN
 !
-!--       Set pavement depth to nzt_soil. Please note, this is just a 
-!--       workaround at the moment. 
+!--       Set pavement depth to nzt_soil. Please note, this is just a
+!--       workaround at the moment.
           DO  m = 1, surf_lsm_h%ns
              IF ( surf_lsm_h%pavement_surface(m) )  THEN
@@ -4125 +4124 @@
                 surf_lsm_h%nzt_pavement(m) = nzt_soil
 
-                DO  k = nzb_soil, nzt_soil  
+                DO  k = nzb_soil, nzt_soil
                    surf_lsm_h%lambda_h_def(k,m) =                              &
                        pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
@@ -4139 +4138 @@
 
                    i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
                    j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                                surf_lsm_v(l)%building_covered(m) ) 
+                                                surf_lsm_v(l)%building_covered(m) )
 
                    surf_lsm_v(l)%nzt_pavement(m) = nzt_soil
 
-                   DO  k = nzb_soil, nzt_soil  
+                   DO  k = nzb_soil, nzt_soil
                       surf_lsm_v(l)%lambda_h_def(k,m) =                        &
                        pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
@@ -4161 +4160 @@
        IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 !
-!--       First, initialize soil temperature and moisture. 
-!--       According to the initialization for surface and soil parameters, 
-!--       initialize soil moisture and temperature via a level approach. This 
-!--       is to assure that all surface elements are initialized, even if 
+!--       First, initialize soil temperature and moisture.
+!--       According to the initialization for surface and soil parameters,
+!--       initialize soil moisture and temperature via a level approach. This
+!--       is to assure that all surface elements are initialized, even if
 !--       data provided from input file contains fill values at some locations.
 !--       Level 1, initialization via profiles given in parameter file
@@ -4170 +4169 @@
              IF ( surf_lsm_h%vegetation_surface(m)  .OR.                       &
                   surf_lsm_h%pavement_surface(m) )  THEN
-                DO  k = nzb_soil, nzt_soil  
+                DO  k = nzb_soil, nzt_soil
                    t_soil_h%var_2d(k,m) = soil_temperature(k)
                    m_soil_h%var_2d(k,m) = soil_moisture(k)
@@ -4181 +4180 @@
                 IF ( surf_lsm_v(l)%vegetation_surface(m)  .OR.                 &
                      surf_lsm_v(l)%pavement_surface(m) )  THEN
-                   DO  k = nzb_soil, nzt_soil  
+                   DO  k = nzb_soil, nzt_soil
                       t_soil_v(l)%var_2d(k,m) = soil_temperature(k)
                       m_soil_v(l)%var_2d(k,m) = soil_moisture(k)
@@ -4190 +4189 @@
           ENDDO
 !
-!--       Level 2 initialization of the soil, read soil properties from 
-!--       dynamic input file. 
+!--       Level 2 initialization of the soil, read soil properties from
+!--       dynamic input file.
           IF ( input_pids_dynamic )  THEN
 !
@@ -4208 +4207 @@
              ALLOCATE( vars_pids(1:num_var_pids) )
              CALL inquire_variable_names( pids_id, vars_pids )
-!           
-!--          Read vertical dimension for soil depth. 
+!
+!--          Read vertical dimension for soil depth.
              IF ( check_existence( vars_pids, 'zsoil' ) )                      &
                 CALL get_dimension_length( pids_id, init_3d%nzs, 'zsoil' )
-!           
+!
 !--          Read also the horizontal dimensions required for soil initialization.
-!--          Please note, in case of non-nested runs or in case of root domain, 
+!--          Please note, in case of non-nested runs or in case of root domain,
 !--          these data is already available, but will be read again for the sake
-!--          of clearness. 
+!--          of clearness.
              CALL get_dimension_length( pids_id, init_3d%nx, 'x'  )
              CALL get_dimension_length( pids_id, init_3d%ny, 'y'  )
-!           
+!
 !--          Check for correct horizontal and vertical dimension. Please note,
-!--          in case of non-nested runs or in case of root domain, these checks 
+!--          in case of non-nested runs or in case of root domain, these checks
 !--          are already performed
              IF ( init_3d%nx-1 /= nx  .OR.  init_3d%ny-1 /= ny )  THEN
@@ -4229 +4228 @@
                 CALL message( 'lsm_init', 'PA0543', 1, 2, 0, 6, 0 )
              ENDIF
-!           
+!
 !--          Read vertical dimensions. Later, these are required for eventual
 !--          inter- and extrapolations of the initialization data.
@@ -4236 +4235 @@
                 CALL get_variable( pids_id, 'zsoil', init_3d%z_soil )
              ENDIF
-!           
+!
 !--          Read initial data for soil moisture
              IF ( check_existence( vars_pids, 'init_soil_m' ) )  THEN
-!           
+!
 !--             Read attributes for the fill value and level-of-detail
                 CALL get_attribute( pids_id, char_fill,                        &
@@ -4247 +4246 @@
                                     init_3d%lod_msoil,                         &
                                     .FALSE., 'init_soil_m' )
-!           
+!
 !--             level-of-detail 1 - read initialization profile
                 IF ( init_3d%lod_msoil == 1 )  THEN
                    ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
-            
+
                    CALL get_variable( pids_id, 'init_soil_m',                  &
                                       init_3d%msoil_1d(0:init_3d%nzs-1) )
-!           
+!
 !--             level-of-detail 2 - read 3D initialization data
                 ELSEIF ( init_3d%lod_msoil == 2 )  THEN
                    ALLOCATE ( init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
-            
-                  CALL get_variable( pids_id, 'init_soil_m',                   &   
+
+                  CALL get_variable( pids_id, 'init_soil_m',                   &
                              init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
                              nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
-            
+
                 ENDIF
                 init_3d%from_file_msoil = .TRUE.
              ENDIF
-!           
+!
 !--          Read soil temperature
              IF ( check_existence( vars_pids, 'init_soil_t' ) )  THEN
-!           
+!
 !--             Read attributes for the fill value and level-of-detail
                 CALL get_attribute( pids_id, char_fill,                        &
@@ -4277 +4276 @@
                                     init_3d%lod_tsoil,                         &
                                     .FALSE., 'init_soil_t' )
-!           
+!
 !--             level-of-detail 1 - read initialization profile
                 IF ( init_3d%lod_tsoil == 1 )  THEN
                    ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
-            
+
                    CALL get_variable( pids_id, 'init_soil_t',                  &
                                       init_3d%tsoil_1d(0:init_3d%nzs-1) )
-            
-!           
+
+!
 !--             level-of-detail 2 - read 3D initialization data
                 ELSEIF ( init_3d%lod_tsoil == 2 )  THEN
                    ALLOCATE ( init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
-                   
-                   CALL get_variable( pids_id, 'init_soil_t',                  &   
+
+                   CALL get_variable( pids_id, 'init_soil_t',                  &
                              init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
                              nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
@@ -4296 +4295 @@
                 init_3d%from_file_tsoil = .TRUE.
              ENDIF
-!           
+!
 !--          Close the input file and deallocate temporary arrays
              DEALLOCATE( vars_pids )
 
              CALL close_input_file( pids_id )
-#endif      
-!           
+#endif
+!
 !--          End of CPU measurement
              CALL cpu_log( log_point_s(85), 'NetCDF input init', 'stop' )
           ENDIF
 !
-!--       In case no dynamic input is available for a child domain but the 
-!--       parent is initialized with dynamic input file, the different soil 
-!--       states can lead to significant discrepancies in the atmospheric 
+!--       In case no dynamic input is available for a child domain but the
+!--       parent is initialized with dynamic input file, the different soil
+!--       states can lead to significant discrepancies in the atmospheric
 !--       surface forcing. For this reason, the child domain is initialized with
 !--       domain-averaged soil profiles from the root domain, even if
-!--       no initialization with inifor is set. Note, as long as a dynamic 
+!--       no initialization with inifor is set. Note, as long as a dynamic
 !--       input file with soil information is available for the child domain,
 !--       the input file information will be used.
@@ -4318 +4317 @@
 #if defined( __parallel )
 !
-!--          Check if soil moisture and temperature in the root model are 
-!--          initialized from dynamic input. This case, distribute these 
+!--          Check if soil moisture and temperature in the root model are
+!--          initialized from dynamic input. This case, distribute these
 !--          information to its child domain(s).
              IF ( pmc_is_rootmodel() )  THEN
-                init_msoil_from_driver_root = init_3d%from_file_msoil 
+                init_msoil_from_driver_root = init_3d%from_file_msoil
                 init_tsoil_from_driver_root = init_3d%from_file_tsoil
              ENDIF
@@ -4331 +4330 @@
                              0, MPI_COMM_WORLD, ierr )
 !
-!--          In case of a nested run, first average the soil profiles in the 
+!--          In case of a nested run, first average the soil profiles in the
 !--          root domain.
              IF ( init_msoil_from_driver_root  .OR.                            &
@@ -4338 +4337 @@
                 IF ( pmc_is_rootmodel() )  THEN
 !
-!--                Child domains will be only initialized with horizontally 
-!--                averaged soil profiles in parent domain (for sake of simplicity). 
-!--                If required, average soil data on root parent domain before the 
+!--                Child domains will be only initialized with horizontally
+!--                averaged soil profiles in parent domain (for sake of simplicity).
+!--                If required, average soil data on root parent domain before the
 !--                soil profiles are distributed onto the child domains.
 !--                Start with soil moisture.
@@ -4349 +4348 @@
                          pr_soil_init(k) = SUM( init_3d%msoil_3d(k,nys:nyn,nxl:nxr)  )
                       ENDDO
-!               
-!--                   Allocate 1D array for soil-moisture profile (will not be 
-!--                   allocated in lod==2 case). 
+!
+!--                   Allocate 1D array for soil-moisture profile (will not be
+!--                   allocated in lod==2 case).
                       ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
                       init_3d%msoil_1d = 0.0_wp
@@ -4357 +4356 @@
                                           SIZE(pr_soil_init),                  &
                                           MPI_REAL, MPI_SUM, comm2d, ierr )
-                
+
                       init_3d%msoil_1d = init_3d%msoil_1d /                    &
                                         REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
                       DEALLOCATE( pr_soil_init )
                    ENDIF
-!               
+!
 !--                Proceed with soil temperature.
                    IF ( init_3d%from_file_tsoil  .AND.                         &
                         init_3d%lod_tsoil == 2 )  THEN
                       ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
-                
+
                       DO  k = 0, init_3d%nzs-1
                          pr_soil_init(k) = SUM( init_3d%tsoil_3d(k,nys:nyn,nxl:nxr)  )
                       ENDDO
-!               
-!--                   Allocate 1D array for soil-temperature profile (will not be 
-!--                   allocated in lod==2 case). 
+!
+!--                   Allocate 1D array for soil-temperature profile (will not be
+!--                   allocated in lod==2 case).
                       ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
                       init_3d%tsoil_1d = 0.0_wp
@@ -4382 +4381 @@
                                         REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
                       DEALLOCATE( pr_soil_init )
-                
+
                    ENDIF
                 ENDIF
-!               
-!--             Broadcast number of soil layers in root model to all childs. 
-!--             Note, only process 0 in COMM_WORLD is sending. 
+!
+!--             Broadcast number of soil layers in root model to all childs.
+!--             Note, only process 0 in COMM_WORLD is sending.
                 IF ( pmc_is_rootmodel() )  nzs_root = init_3d%nzs
-                
+
                 CALL MPI_BCAST( nzs_root, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
-!               
-!--             Allocate dummy arrays for soil moisture and temperature profiles 
-!--             on all domains.             
+!
+!--             Allocate dummy arrays for soil moisture and temperature profiles
+!--             on all domains.
                 ALLOCATE( z_soil_root(1:nzs_root)   )
                 IF ( init_msoil_from_driver_root )                             &
@@ -4421 +4420 @@
 !
 !--             In the following, the child domains decide whether they take
-!--             the information from the root domain or not. 
+!--             the information from the root domain or not.
                 IF ( .NOT. pmc_is_rootmodel() )  THEN
-!               
+!
 !--                If soil moisture or temperature isn't in dynamic input file for
-!--                the child, take the information provided from the root model. 
+!--                the child, take the information provided from the root model.
 !--                Start with z-dimension
                    IF ( .NOT. init_3d%from_file_msoil  .OR.                    &
@@ -4433 +4432 @@
                       init_3d%z_soil(1:init_3d%nzs) = z_soil_root
                    ENDIF
-!                  
-!--                Take soil moisture. Note, control flags from_file... and LoD 
-!--                need to be set. 
+!
+!--                Take soil moisture. Note, control flags from_file... and LoD
+!--                need to be set.
                    IF ( .NOT. init_3d%from_file_msoil )  THEN
                       ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
                       init_3d%lod_msoil = 1
                       init_3d%from_file_msoil = .TRUE.
-                      
-                      init_3d%msoil_1d = m_soil_root             
+
+                      init_3d%msoil_1d = m_soil_root
                    ENDIF
-!                  
-!--                Take soil temperature. Note, control flags from_file... and LoD 
-!--                need to be set. 
+!
+!--                Take soil temperature. Note, control flags from_file... and LoD
+!--                need to be set.
                    IF (  .NOT. init_3d%from_file_tsoil )  THEN
                       ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
                       init_3d%lod_tsoil = 1
                       init_3d%from_file_tsoil = .TRUE.
-                      
-                      init_3d%tsoil_1d = t_soil_root  
+
+                      init_3d%tsoil_1d = t_soil_root
                    ENDIF
                 ENDIF
-                
+
                 DEALLOCATE( z_soil_root )
                 DEALLOCATE( m_soil_root )
@@ -4462 +4461 @@
           ENDIF
 !
-!--       Proceed with Level 2 initialization. 
+!--       Proceed with Level 2 initialization.
           IF ( init_3d%from_file_msoil )  THEN
 
@@ -4516 +4515 @@
 !--                      input do not need to be checked whether a grid point
 !--                      is building covered. This is because soil data in the
-!--                      dynamic input is provided for the whole domain.  
+!--                      dynamic input is provided for the whole domain.
                          i = surf_lsm_v(l)%i(m)
                          j = surf_lsm_v(l)%j(m)
-                         
+
                          IF ( init_3d%msoil_3d(0,j,i) /= init_3d%fill_msoil )  &
                             CALL interpolate_soil_profile(                     &
@@ -4575 +4574 @@
                       i = surf_lsm_h%i(m)
                       j = surf_lsm_h%j(m)
-                      
+
                       IF ( init_3d%tsoil_3d(0,j,i) /= init_3d%fill_tsoil )     &
                          CALL interpolate_soil_profile(                        &
@@ -4596 +4595 @@
 !--                      input do not need to be checked whether a grid point
 !--                      is building covered. This is because soil data in the
-!--                      dynamic input is provided for the whole domain.  
+!--                      dynamic input is provided for the whole domain.
                          i = surf_lsm_v(l)%i(m)
                          j = surf_lsm_v(l)%j(m)
-                         
+
                          IF ( init_3d%tsoil_3d(0,j,i) /= init_3d%fill_tsoil )  &
                             CALL interpolate_soil_profile(                     &
@@ -4617 +4616 @@
           ENDIF
 !
-!--       After soil moisture and temperature are finally initialized, check 
-!--       if soil moisture is higher than its saturation value. If this would 
-!--       be the case, the soil model parametrization will produce floating 
+!--       After soil moisture and temperature are finally initialized, check
+!--       if soil moisture is higher than its saturation value. If this would
+!--       be the case, the soil model parametrization will produce floating
 !--       point errors. Hence, limit the soil moisture to its saturation value
-!--       and give a warning. 
+!--       and give a warning.
           DO  m = 1, surf_lsm_h%ns
              IF ( surf_lsm_h%vegetation_surface(m)  .OR.                       &
@@ -4633 +4632 @@
                             'thus limited to this value to maintain stability.'
                       CALL message( 'lsm_init', 'PA0458', 0, 1, myid, 6, 0 )
-                   ENDIF               
+                   ENDIF
                 ENDDO
              ENDIF
@@ -4651 +4650 @@
                             ' and is ' //                                      &
                             'thus limited to this value to maintain stability.'
-                         CALL message( 'lsm_init', 'PA0458', 0, 1, myid, 6, 0 )                   
+                         CALL message( 'lsm_init', 'PA0458', 0, 1, myid, 6, 0 )
                       ENDIF
                    ENDDO
@@ -4662 +4661 @@
           DO  m = 1, surf_lsm_h%ns
 
-             i   = surf_lsm_h%i(m)            
+             i   = surf_lsm_h%i(m)
              j   = surf_lsm_h%j(m)
              k   = surf_lsm_h%k(m)
@@ -4670 +4669 @@
              t_surface_h%var_1d(m)    = t_soil_h%var_2d(nzb_soil,m)
              surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exner(nzb)
-             
+
              IF ( bulk_cloud_model  .OR. cloud_droplets ) THEN
                 surf_lsm_h%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
              ELSE
                 surf_lsm_h%pt1(m) = pt(k,j,i)
-             ENDIF 
+             ENDIF
 !
 !--          Assure that r_a cannot be zero at model start
@@ -4691 +4690 @@
           DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
-                i   = surf_lsm_v(l)%i(m)            
+                i   = surf_lsm_v(l)%i(m)
                 j   = surf_lsm_v(l)%j(m)
-                k   = surf_lsm_v(l)%k(m)         
+                k   = surf_lsm_v(l)%k(m)
 !
 !--             Initialize surface temperature with soil temperature in the uppermost
@@ -4704 +4703 @@
                 ELSE
                    surf_lsm_v(l)%pt1(m) = pt(k,j,i)
-                ENDIF 
+                ENDIF
 
 !
@@ -4711 +4710 @@
                      surf_lsm_v(l)%pt1(m) = surf_lsm_v(l)%pt1(m) + 1.0E-20_wp
 !
-!--             Set artifical values for ts and us so that r_a has its initial value 
+!--             Set artifical values for ts and us so that r_a has its initial value
 !--             for the first time step. Only for interior core domain, not for ghost points
                 surf_lsm_v(l)%us(m)   = 0.1_wp
@@ -4798 +4797 @@
 !--       Map calculated root fractions
           DO  m = 1, surf_lsm_h%ns
-             DO  k = nzb_soil, nzt_soil 
+             DO  k = nzb_soil, nzt_soil
                 IF ( surf_lsm_h%pavement_surface(m)  .AND.                     &
                      k <= surf_lsm_h%nzt_pavement(m) )  THEN
@@ -4808 +4807 @@
              ENDDO
 !
-!--          Normalize so that the sum = 1. Only relevant when the root         
+!--          Normalize so that the sum = 1. Only relevant when the root
 !--          distribution was set to zero due to pavement at some layers.
              IF ( SUM( surf_lsm_h%root_fr(:,m) ) > 0.0_wp )  THEN
@@ -4828 +4827 @@
                 ENDDO
 !
-!--             Normalize so that the sum = 1. Only relevant when the root      
+!--             Normalize so that the sum = 1. Only relevant when the root
 !--             distribution was set to zero due to pavement at some layers.
                 IF ( SUM( surf_lsm_v(l)%root_fr(:,m) ) > 0.0_wp )  THEN
-                   DO  k = nzb_soil, nzt_soil  
+                   DO  k = nzb_soil, nzt_soil
                       surf_lsm_v(l)%root_fr(k,m) = surf_lsm_v(l)%root_fr(k,m)  &
                       / SUM( surf_lsm_v(l)%root_fr(:,m) )
@@ -4846 +4845 @@
              IF ( surf_lsm_h%vegetation_surface(m) )  THEN
                 i = surf_lsm_h%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,            &
-                                             surf_lsm_v(l)%building_covered(m) ) 
+                                             surf_lsm_v(l)%building_covered(m) )
                 j = surf_lsm_h%j(m) + MERGE( 0, surf_lsm_v(l)%joff,            &
-                                             surf_lsm_v(l)%building_covered(m) ) 
-                DO  k = nzb_soil, nzt_soil 
-                   surf_lsm_h%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i) 
+                                             surf_lsm_v(l)%building_covered(m) )
+                DO  k = nzb_soil, nzt_soil
+                   surf_lsm_h%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i)
                 ENDDO
 
@@ -4860 +4859 @@
                 IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
                    i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
+                                                   surf_lsm_v(l)%building_covered(m) )
                    j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                                   surf_lsm_v(l)%building_covered(m) ) 
-
-                   DO  k = nzb_soil, nzt_soil 
-                      surf_lsm_v(l)%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i) 
+                                                   surf_lsm_v(l)%building_covered(m) )
+
+                   DO  k = nzb_soil, nzt_soil
+                      surf_lsm_v(l)%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i)
                    ENDDO
 
@@ -4873 +4872 @@
 
        ENDIF
- 
+
 !
 !--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
@@ -4880 +4879 @@
 
 !
-!--    Calculate new roughness lengths (for water surfaces only, i.e. only 
+!--    Calculate new roughness lengths (for water surfaces only, i.e. only
 !-     horizontal surfaces)
        IF ( .NOT. constant_roughness )  CALL calc_z0_water_surface
@@ -4896 +4895 @@
 
 
-!--    Store initial profiles of t_soil and m_soil (assuming they are 
+!--    Store initial profiles of t_soil and m_soil (assuming they are
 !--    horizontally homogeneous on this PE)
 !--    DEACTIVATED FOR NOW - leads to error when number of locations with
@@ -4902 +4901 @@
 !        hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
 !                                                 2, statistic_regions+1 )
-!        hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1),  & 
+!        hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
 !                                                 2, statistic_regions+1 )
 
 !
-!--    Finally, make some consistency checks. 
+!--    Finally, make some consistency checks.
 !--    Ceck for illegal combination of LAI and vegetation coverage.
        IF ( ANY( .NOT. surf_lsm_h%pavement_surface  .AND.                      &
@@ -4926 +4925 @@
        ENDDO
 !
-!--    Check if roughness length for momentum, heat, or moisture exceed 
-!--    surface-layer height and decrease local roughness length where 
+!--    Check if roughness length for momentum, heat, or moisture exceed
+!--    surface-layer height and decrease local roughness length where
 !--    necessary. This case, give an informative message. Note, to avoid
 !--    that the job-protocoll is messed-up, this message is only given once.
@@ -5011 +5010 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE lsm_init_arrays
-    
+
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  l !< index indicating facing of surface array 
-   
+       INTEGER(iwp) ::  l !< index indicating facing of surface array
+
        ALLOCATE ( root_extr(nzb_soil:nzt_soil) )
-       root_extr = 0.0_wp 
-       
-!
-!--    Allocate surface and soil temperature / humidity. Please note, 
+       root_extr = 0.0_wp
+
+!
+!--    Allocate surface and soil temperature / humidity. Please note,
 !--    these arrays are allocated according to surface-data structure,
-!--    even if they do not belong to the data type due to the 
-!--    pointer arithmetric (TARGET attribute is not allowed in a data-type). 
+!--    even if they do not belong to the data type due to the
+!--    pointer arithmetric (TARGET attribute is not allowed in a data-type).
 !
 !--    Horizontal surfaces
@@ -5063 +5062 @@
        ALLOCATE ( tt_surface_h_m%var_1d(1:surf_lsm_h%ns)                 )
        ALLOCATE ( tm_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
-       ALLOCATE ( tt_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  ) 
+       ALLOCATE ( tt_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
 !
 !--    Horizontal surfaces
@@ -5071 +5070 @@
           ALLOCATE ( tm_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
           ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
-       ENDDO 
+       ENDDO
 
 !
@@ -5091 +5090 @@
        ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns)            )
        ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns)            )
-       ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns)           ) 
+       ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns)           )
        ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns)                 )
        ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)            )
@@ -5103 +5102 @@
        surf_lsm_h%water_surface        = .FALSE.
        surf_lsm_h%pavement_surface     = .FALSE.
-       surf_lsm_h%vegetation_surface   = .FALSE. 
+       surf_lsm_h%vegetation_surface   = .FALSE.
 
 !
@@ -5139 +5138 @@
           surf_lsm_v(l)%water_surface       = .FALSE.
           surf_lsm_v(l)%pavement_surface    = .FALSE.
-          surf_lsm_v(l)%vegetation_surface  = .FALSE. 
-          
+          surf_lsm_v(l)%vegetation_surface  = .FALSE.
+
 
 !
 !--       Set default values
           surf_lsm_v(l)%r_canopy_min = 0.0_wp
-        
+
        ENDDO
 
@@ -5178 +5177 @@
        IMPLICIT NONE
 
-       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
-       
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
+
        NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
                                   canopy_resistance_coefficient,               &
@@ -5186 +5185 @@
                                   deep_soil_temperature,                       &
                                   dz_soil,                                     &
-                                  f_shortwave_incoming, field_capacity,        & 
+                                  f_shortwave_incoming, field_capacity,        &
                                   aero_resist_kray, hydraulic_conductivity,    &
                                   lambda_surface_stable,                       &
@@ -5214 +5213 @@
                                   deep_soil_temperature,                       &
                                   dz_soil,                                     &
-                                  f_shortwave_incoming, field_capacity,        & 
+                                  f_shortwave_incoming, field_capacity,        &
                                   aero_resist_kray, hydraulic_conductivity,    &
                                   lambda_surface_stable,                       &
@@ -5234 +5233 @@
                                   z0h_water, z0q_water, z0_pavement,           &
                                   z0h_pavement, z0q_pavement
-                                  
+
        line = ' '
-  
+
 !
 !--    Try to find land surface model package
@@ -5253 +5252 @@
 !--    Set flag that indicates that the land surface model is switched on
        land_surface = .TRUE.
-       
+
        GOTO 14
 
@@ -5276 +5275 @@
                      'land_surface_parameters instead'
        CALL message( 'lsm_parin', 'PA0487', 0, 1, 0, 6, 0 )
-       
+
 !
 !--    Set flag that indicates that the land surface model is switched on
        land_surface = .TRUE.
-       
+
        GOTO 14
 
@@ -5289 +5288 @@
 
  14    CONTINUE
-       
+
 
     END SUBROUTINE lsm_parin
@@ -5363 +5362 @@
                 IF ( surf%pavement_surface(m)  .AND.                           &
                      k <= surf%nzt_pavement(m) )  THEN
-                   
+
                    surf%rho_c_total(k,m) = surf%rho_c_total_def(k,m)
-                   lambda_temp(k)        = surf%lambda_h_def(k,m) 
-
-                ELSE            
-!
-!--                Calculate volumetric heat capacity of the soil, taking 
-!--                into account water content 
+                   lambda_temp(k)        = surf%lambda_h_def(k,m)
+
+                ELSE
+!
+!--                Calculate volumetric heat capacity of the soil, taking
+!--                into account water content
                    surf%rho_c_total(k,m) = (rho_c_soil *                       &
                                                ( 1.0_wp - surf%m_sat(k,m) )    &
@@ -5390 +5389 @@
 
 !
-!--          Calculate soil heat conductivity (lambda_h) at the _layer level 
+!--          Calculate soil heat conductivity (lambda_h) at the _layer level
 !--          using linear interpolation. For pavement surface, the
 !--          true pavement depth is considered
@@ -5423 +5422 @@
                                        surf_t_soil%var_2d(nzb_soil:nzt_soil,m) &
                                                + dt_3d * ( tsc(2)              &
-                                               * tend(nzb_soil:nzt_soil)       & 
+                                               * tend(nzb_soil:nzt_soil)       &
                                                + tsc(3)                        &
                                                * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
@@ -5448 +5447 @@
 
 !
-!--             In order to prevent water tranport through paved surfaces, 
+!--             In order to prevent water tranport through paved surfaces,
 !--             conductivity and diffusivity are set to zero
                 IF ( surf%pavement_surface(m)  .AND.                           &
@@ -5454 +5453 @@
                    lambda_temp(k) = 0.0_wp
                    gamma_temp(k)  = 0.0_wp
-   
-                ELSE 
-    
+
+                ELSE
+
 !
 !--                Calculate soil diffusivity at the center of the soil layers
@@ -5500 +5499 @@
                 IF ( humidity )  THEN
 !
-!--                Calculate soil diffusivity (lambda_w) at the _layer level 
+!--                Calculate soil diffusivity (lambda_w) at the _layer level
 !--                using linear interpolation. To do: replace this with
 !--                ECMWF-IFS Eq. 8.81
                    DO  k = nzb_soil, nzt_soil-1
-                
+
                       surf%lambda_w(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
                                            * 0.5_wp
                       surf%gamma_w(k,m)  = ( gamma_temp(k+1)  +  gamma_temp(k) )  &
                                            * 0.5_wp
-                                            
+
                    ENDDO
 !
 !
-!--                In case of a closed bottom (= water content is conserved), 
-!--                set hydraulic conductivity to zero to that no water will be 
+!--                In case of a closed bottom (= water content is conserved),
+!--                set hydraulic conductivity to zero to that no water will be
 !--                lost in the bottom layer. As gamma_w is always a positive value,
 !--                it cannot be set to zero in case of purely dry soil since this
@@ -5525 +5524 @@
                    ELSE
                       surf%gamma_w(nzt_soil,m) = gamma_temp(nzt_soil)
-                   ENDIF     
-
-!--                The root extraction (= root_extr * qsws_veg / (rho_l     
-!--                * l_v)) ensures the mass conservation for water. The         
-!--                transpiration of plants equals the cumulative withdrawals by 
-!--                the roots in the soil. The scheme takes into account the 
-!--                availability of water in the soil layers as well as the root 
-!--                fraction in the respective layer. Layer with moisture below 
-!--                wilting point will not contribute, which reflects the 
+                   ENDIF
+
+!--                The root extraction (= root_extr * qsws_veg / (rho_l
+!--                * l_v)) ensures the mass conservation for water. The
+!--                transpiration of plants equals the cumulative withdrawals by
+!--                the roots in the soil. The scheme takes into account the
+!--                availability of water in the soil layers as well as the root
+!--                fraction in the respective layer. Layer with moisture below
+!--                wilting point will not contribute, which reflects the
 !--                preference of plants to take water from moister layers.
 !
-!--                Calculate the root extraction (ECMWF 7.69, the sum of 
-!--                root_extr = 1). The energy balance solver guarantees a 
-!--                positive transpiration, so that there is no need for an 
+!--                Calculate the root extraction (ECMWF 7.69, the sum of
+!--                root_extr = 1). The energy balance solver guarantees a
+!--                positive transpiration, so that there is no need for an
 !--                additional check.
                    m_total = 0.0_wp
@@ -5546 +5545 @@
                                 * surf_m_soil%var_2d(k,m)
                       ENDIF
-                   ENDDO  
+                   ENDDO
                    IF ( m_total > 0.0_wp )  THEN
                       DO  k = nzb_soil, nzt_soil
@@ -5587 +5586 @@
                                    root_extr(nzt_soil)                         &
                                    * surf%qsws_veg(m) * drho_l_lv )            &
-                                  ) * ddz_soil(nzt_soil)             
+                                  ) * ddz_soil(nzt_soil)
 
                    surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) =                 &
                                        surf_m_soil%var_2d(nzb_soil:nzt_soil,m) &
                                          + dt_3d * ( tsc(2) * tend(:)          &
-                                         + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
-   
+                                         + tsc(3) * surf_tm_soil_m%var_2d(:,m) )
+
 !
 !--                Account for dry and wet soils to keep solution stable
@@ -5599 +5598 @@
                    DO  k = nzb_soil, nzt_soil
                       surf_m_soil_p%var_2d(k,m) = MIN( surf_m_soil_p%var_2d(k,m), surf%m_sat(k,m) )
-                      surf_m_soil_p%var_2d(k,m) = MAX( surf_m_soil_p%var_2d(k,m), 0.0_wp )                      
+                      surf_m_soil_p%var_2d(k,m) = MAX( surf_m_soil_p%var_2d(k,m), 0.0_wp )
                    ENDDO
-  
+
 !
 !--                Calculate m_soil tendencies for the next Runge-Kutta step
@@ -5618 +5617 @@
 
                       ENDIF
-                      
+
                    ENDIF
-                   
+
                 ENDIF
 
@@ -5638 +5637 @@
     END SUBROUTINE lsm_soil_model
 
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -5650 +5649 @@
        INTEGER, INTENT(IN) :: mod_count
 
-    
+
        SELECT CASE ( mod_count )
 
@@ -5708 +5707 @@
 !------------------------------------------------------------------------------!
 SUBROUTINE lsm_3d_data_averaging( mode, variable )
- 
+
 
     USE control_parameters
@@ -5716 +5715 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  mode    !< 
-    CHARACTER (LEN=*) :: variable !< 
-
-    INTEGER(iwp) ::  i       !< 
-    INTEGER(iwp) ::  j       !< 
-    INTEGER(iwp) ::  k       !< 
+    CHARACTER (LEN=*) ::  mode    !<
+    CHARACTER (LEN=*) :: variable !<
+
+    INTEGER(iwp) ::  i       !<
+    INTEGER(iwp) ::  j       !<
+    INTEGER(iwp) ::  k       !<
     INTEGER(iwp) ::  m       !< running index
 
@@ -5804 +5803 @@
 
           CASE ( 'c_liq*' )
-             IF ( ALLOCATED( c_liq_av ) ) THEN 
+             IF ( ALLOCATED( c_liq_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    c_liq_av(j,i) = c_liq_av(j,i) + surf_lsm_h%c_liq(m)
                 ENDDO
-             ENDIF   
+             ENDIF
 
           CASE ( 'c_soil*' )
-             IF ( ALLOCATED( c_soil_av ) ) THEN 
+             IF ( ALLOCATED( c_soil_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - surf_lsm_h%c_veg(m))
@@ -5822 +5821 @@
 
           CASE ( 'c_veg*' )
-             IF ( ALLOCATED( c_veg_av ) ) THEN 
+             IF ( ALLOCATED( c_veg_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    c_veg_av(j,i) = c_veg_av(j,i) + surf_lsm_h%c_veg(m)
@@ -5831 +5830 @@
 
           CASE ( 'lai*' )
-             IF ( ALLOCATED( lai_av ) ) THEN 
+             IF ( ALLOCATED( lai_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    lai_av(j,i) = lai_av(j,i) + surf_lsm_h%lai(m)
@@ -5840 +5839 @@
 
           CASE ( 'm_liq*' )
-             IF ( ALLOCATED( m_liq_av ) ) THEN 
+             IF ( ALLOCATED( m_liq_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    m_liq_av(j,i) = m_liq_av(j,i) + m_liq_h%var_1d(m)
@@ -5849 +5848 @@
 
           CASE ( 'm_soil' )
-             IF ( ALLOCATED( m_soil_av ) ) THEN 
+             IF ( ALLOCATED( m_soil_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    DO  k = nzb_soil, nzt_soil
@@ -5860 +5859 @@
 
           CASE ( 'qsws_liq*' )
-             IF ( ALLOCATED( qsws_liq_av ) ) THEN 
+             IF ( ALLOCATED( qsws_liq_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    qsws_liq_av(j,i) = qsws_liq_av(j,i) +                       &
@@ -5870 +5869 @@
 
           CASE ( 'qsws_soil*' )
-             IF ( ALLOCATED( qsws_soil_av ) ) THEN 
+             IF ( ALLOCATED( qsws_soil_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    qsws_soil_av(j,i) = qsws_soil_av(j,i) +                     &
@@ -5880 +5879 @@
 
           CASE ( 'qsws_veg*' )
-             IF ( ALLOCATED(qsws_veg_av ) ) THEN 
+             IF ( ALLOCATED(qsws_veg_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    qsws_veg_av(j,i) = qsws_veg_av(j,i) +                       &
@@ -5890 +5889 @@
 
           CASE ( 'r_s*' )
-             IF ( ALLOCATED( r_s_av) ) THEN 
+             IF ( ALLOCATED( r_s_av) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    r_s_av(j,i) = r_s_av(j,i) + surf_lsm_h%r_s(m)
@@ -5899 +5898 @@
 
           CASE ( 't_soil' )
-             IF ( ALLOCATED( t_soil_av ) ) THEN 
+             IF ( ALLOCATED( t_soil_av ) ) THEN
                 DO  m = 1, surf_lsm_h%ns
-                   i   = surf_lsm_h%i(m)            
+                   i   = surf_lsm_h%i(m)
                    j   = surf_lsm_h%j(m)
                    DO  k = nzb_soil, nzt_soil
@@ -5908 +5907 @@
                 ENDDO
              ENDIF
-             
+
           CASE DEFAULT
              CONTINUE
@@ -6014 +6013 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      r_s_av(j,i) = r_s_av(j,i)                                & 
+                      r_s_av(j,i) = r_s_av(j,i)                                &
                                     / REAL( average_count_3d, KIND=wp )
                    ENDDO
@@ -6032 +6031 @@
              ENDIF
 !
-!-- 
+!--
 
        END SELECT
@@ -6049 +6048 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
-    
+
      IMPLICIT NONE
 
-     CHARACTER (LEN=*), INTENT(IN)  ::  var         !< 
-     LOGICAL, INTENT(OUT)           ::  found       !< 
-     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
-     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
-     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
+     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
+     LOGICAL, INTENT(OUT)           ::  found       !<
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
 
      found  = .TRUE.
@@ -6087 +6086 @@
  SUBROUTINE lsm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
                                 two_d, nzb_do, nzt_do )
- 
+
     USE indices
 
@@ -6093 +6092 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  grid     !< 
-    CHARACTER (LEN=*) ::  mode     !< 
-    CHARACTER (LEN=*) ::  variable !< 
-
-    INTEGER(iwp) ::  av      !< 
-    INTEGER(iwp) ::  i       !< running index 
+    CHARACTER (LEN=*) ::  grid     !<
+    CHARACTER (LEN=*) ::  mode     !<
+    CHARACTER (LEN=*) ::  variable !<
+
+    INTEGER(iwp) ::  av      !<
+    INTEGER(iwp) ::  i       !< running index
     INTEGER(iwp) ::  j       !< running index
     INTEGER(iwp) ::  k       !< running index
     INTEGER(iwp) ::  m       !< running index
-    INTEGER(iwp) ::  nzb_do  !< 
-    INTEGER(iwp) ::  nzt_do  !< 
-
-    LOGICAL      ::  found !< 
+    INTEGER(iwp) ::  nzb_do  !<
+    INTEGER(iwp) ::  nzt_do  !<
+
+    LOGICAL      ::  found !<
     LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
 
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
 
 
@@ -6117 +6116 @@
     SELECT CASE ( TRIM( variable ) )
 !
-!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein. 
-!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged. 
+!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
+!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
        CASE ( 'c_liq*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) = surf_lsm_h%c_liq(m) * surf_lsm_h%c_veg(m)
@@ -6144 +6143 @@
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) = 1.0_wp - surf_lsm_h%c_veg(m)
@@ -6166 +6165 @@
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) = surf_lsm_h%c_veg(m)
@@ -6188 +6187 @@
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) = surf_lsm_h%lai(m)
@@ -6210 +6209 @@
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) = m_liq_h%var_1d(m)
@@ -6232 +6231 @@
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i   = surf_lsm_h%i(m)            
+                i   = surf_lsm_h%i(m)
                 j   = surf_lsm_h%j(m)
                 DO k = nzb_soil, nzt_soil
@@ -6256 +6255 @@
 
           IF ( mode == 'xy' ) grid = 'zs'
-         
+
        CASE ( 'qsws_liq*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) = surf_lsm_h%qsws_liq(m)
@@ -6270 +6269 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) =  qsws_liq_av(j,i)
                 ENDDO
@@ -6282 +6281 @@
           IF ( av == 0 ) THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_soil(m)
@@ -6292 +6291 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) =  qsws_soil_av(j,i)
                 ENDDO
@@ -6304 +6303 @@
           IF ( av == 0 ) THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_veg(m)
@@ -6314 +6313 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) =  qsws_veg_av(j,i)
                 ENDDO
@@ -6327 +6326 @@
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
+                i                   = surf_lsm_h%i(m)
                 j                   = surf_lsm_h%j(m)
                 local_pf(i,j,nzb+1) = surf_lsm_h%r_s(m)
@@ -6349 +6348 @@
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
-                i   = surf_lsm_h%i(m)            
+                i   = surf_lsm_h%i(m)
                 j   = surf_lsm_h%j(m)
                 DO k = nzb_soil, nzt_soil
@@ -6380 +6379 @@
 
     END SELECT
- 
+
  END SUBROUTINE lsm_data_output_2d
 
@@ -6391 +6390 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_data_output_3d( av, variable, found, local_pf )
- 
+
 
     USE indices
@@ -6398 +6397 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable !< 
-
-    INTEGER(iwp) ::  av    !< 
-    INTEGER(iwp) ::  i     !< 
-    INTEGER(iwp) ::  j     !< 
-    INTEGER(iwp) ::  k     !< 
+    CHARACTER (LEN=*) ::  variable !<
+
+    INTEGER(iwp) ::  av    !<
+    INTEGER(iwp) ::  i     !<
+    INTEGER(iwp) ::  j     !<
+    INTEGER(iwp) ::  k     !<
     INTEGER(iwp) ::  m     !< running index
 
-    LOGICAL      ::  found !< 
+    LOGICAL      ::  found !<
 
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_soil:nzt_soil) ::  local_pf !< 
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_soil:nzt_soil) ::  local_pf !<
 
 
@@ -6424 +6423 @@
          IF ( av == 0 )  THEN
             DO  m = 1, surf_lsm_h%ns
-                i   = surf_lsm_h%i(m)            
+                i   = surf_lsm_h%i(m)
                 j   = surf_lsm_h%j(m)
                 DO  k = nzb_soil, nzt_soil
@@ -6448 +6447 @@
          IF ( av == 0 )  THEN
             DO  m = 1, surf_lsm_h%ns
-               i   = surf_lsm_h%i(m)            
+               i   = surf_lsm_h%i(m)
                j   = surf_lsm_h%j(m)
                DO  k = nzb_soil, nzt_soil
@@ -6482 +6481 @@
 ! Description:
 ! ------------
-!> Write restart data for land surface model. It is necessary to write 
+!> Write restart data for land surface model. It is necessary to write
 !> start_index and end_index several times.
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_wrd_local
-        
+
 
     IMPLICIT NONE
@@ -6544 +6543 @@
        WRITE ( 14 )  qsws_veg_av
     ENDIF
-    
+
     IF ( ALLOCATED( t_soil_av ) )  THEN
        CALL wrd_write_string( 't_soil_av' )
@@ -6558 +6557 @@
     CALL wrd_write_string( 't_soil_h' )
     WRITE ( 14 )  t_soil_h%var_2d
-       
-
-       
+
+
+
     DO  l = 0, 3
 
@@ -6569 +6568 @@
        WRITE ( 14 )  surf_lsm_v(l)%end_index
 
-       WRITE( dum, '(I1)')  l   
+       WRITE( dum, '(I1)')  l
 
        CALL wrd_write_string( 't_soil_v(' // dum // ')' )
        WRITE ( 14 )  t_soil_v(l)%var_2d
-             
+
     ENDDO
 
@@ -6593 +6592 @@
        WRITE ( 14 )  surf_lsm_v(l)%end_index
 
-       WRITE( dum, '(I1)')  l   
+       WRITE( dum, '(I1)')  l
 
        CALL wrd_write_string( 'm_soil_v(' // dum // ')' )
-       WRITE ( 14 )  m_soil_v(l)%var_2d 
-      
+       WRITE ( 14 )  m_soil_v(l)%var_2d
+
     ENDDO
 
@@ -6608 +6607 @@
     CALL wrd_write_string( 'm_liq_h' )
     WRITE ( 14 )  m_liq_h%var_1d
-       
+
     DO  l = 0, 3
 
@@ -6617 +6616 @@
        WRITE ( 14 )  surf_lsm_v(l)%end_index
 
-       WRITE( dum, '(I1)')  l   
+       WRITE( dum, '(I1)')  l
 
        CALL wrd_write_string( 'm_liq_v(' // dum // ')' )
-       WRITE ( 14 )  m_liq_v(l)%var_1d     
-                
+       WRITE ( 14 )  m_liq_v(l)%var_1d
+
     ENDDO
 
@@ -6641 +6640 @@
        WRITE ( 14 )  surf_lsm_v(l)%end_index
 
-       WRITE( dum, '(I1)')  l   
+       WRITE( dum, '(I1)')  l
 
        CALL wrd_write_string( 't_surface_v(' // dum // ')' )
-       WRITE ( 14 )  t_surface_v(l)%var_1d     
-       
+       WRITE ( 14 )  t_surface_v(l)%var_1d
+
     ENDDO
 
@@ -6661 +6660 @@
                           nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,    &
                           nys_on_file, tmp_2d, found )
- 
+
 
     USE control_parameters
-        
+
     USE indices
-    
+
     USE pegrid
 
@@ -6672 +6671 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  k                 !< 
+    INTEGER(iwp) ::  k                 !<
     INTEGER(iwp) ::  l                 !< running index surface orientation
     INTEGER(iwp) ::  ns_h_on_file_lsm  !< number of horizontal surface elements (natural type) on file
-    INTEGER(iwp) ::  nxlc              !< 
-    INTEGER(iwp) ::  nxlf              !< 
+    INTEGER(iwp) ::  nxlc              !<
+    INTEGER(iwp) ::  nxlf              !<
     INTEGER(iwp) ::  nxl_on_file       !< index of left boundary on former local domain
-    INTEGER(iwp) ::  nxrc              !< 
-    INTEGER(iwp) ::  nxrf              !< 
+    INTEGER(iwp) ::  nxrc              !<
+    INTEGER(iwp) ::  nxrf              !<
     INTEGER(iwp) ::  nxr_on_file       !< index of right boundary on former local domain
-    INTEGER(iwp) ::  nync              !< 
-    INTEGER(iwp) ::  nynf              !< 
+    INTEGER(iwp) ::  nync              !<
+    INTEGER(iwp) ::  nynf              !<
     INTEGER(iwp) ::  nyn_on_file       !< index of north boundary on former local domain
-    INTEGER(iwp) ::  nysc              !< 
-    INTEGER(iwp) ::  nysf              !< 
+    INTEGER(iwp) ::  nysc              !<
+    INTEGER(iwp) ::  nysf              !<
     INTEGER(iwp) ::  nys_on_file       !< index of south boundary on former local domain
 
     INTEGER(iwp) ::  ns_v_on_file_lsm(0:3) !< number of vertical surface elements (natural type) on file
 
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  start_index_on_file 
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  start_index_on_file
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  end_index_on_file
 
     LOGICAL, INTENT(OUT)  :: found
 
-    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !< 
+    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
 
     REAL(wp), DIMENSION(nzb_soil:nzt_soil,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
@@ -6701 +6700 @@
     TYPE(surf_type_lsm), SAVE :: tmp_walltype_h_1d   !< temporary 1D array containing the respective surface variable stored on file, horizontal surfaces
     TYPE(surf_type_lsm), SAVE :: tmp_walltype_h_2d   !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces
-    TYPE(surf_type_lsm), SAVE :: tmp_walltype_h_2d2  !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces 
-
-    TYPE(surf_type_lsm), DIMENSION(0:3), SAVE :: tmp_walltype_v_1d   !< temporary 1D array containing the respective surface variable stored on file, vertical surfaces 
-    TYPE(surf_type_lsm), DIMENSION(0:3), SAVE :: tmp_walltype_v_2d   !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces 
+    TYPE(surf_type_lsm), SAVE :: tmp_walltype_h_2d2  !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces
+
+    TYPE(surf_type_lsm), DIMENSION(0:3), SAVE :: tmp_walltype_v_1d   !< temporary 1D array containing the respective surface variable stored on file, vertical surfaces
+    TYPE(surf_type_lsm), DIMENSION(0:3), SAVE :: tmp_walltype_v_2d   !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces
     TYPE(surf_type_lsm), DIMENSION(0:3), SAVE :: tmp_walltype_v_2d2  !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces
 
@@ -6714 +6713 @@
 
        CASE ( 'ns_h_on_file_lsm' )
-          IF ( k == 1 )  THEN  
+          IF ( k == 1 )  THEN
              READ ( 13 ) ns_h_on_file_lsm
 
              IF ( ALLOCATED( tmp_walltype_h_1d%var_1d ) )                      &
-                DEALLOCATE( tmp_walltype_h_1d%var_1d )                         
+                DEALLOCATE( tmp_walltype_h_1d%var_1d )
              IF ( ALLOCATED( tmp_walltype_h_2d%var_2d ) )                      &
-                DEALLOCATE( tmp_walltype_h_2d%var_2d )                         
+                DEALLOCATE( tmp_walltype_h_2d%var_2d )
              IF ( ALLOCATED( tmp_walltype_h_2d2%var_2d ) )                     &
-                DEALLOCATE( tmp_walltype_h_2d2%var_2d )                        
-                                                                               
-!                                                                              
-!--          Allocate temporary arrays to store surface data                   
-             ALLOCATE( tmp_walltype_h_1d%var_1d(1:ns_h_on_file_lsm) )          
+                DEALLOCATE( tmp_walltype_h_2d2%var_2d )
+
+!
+!--          Allocate temporary arrays to store surface data
+             ALLOCATE( tmp_walltype_h_1d%var_1d(1:ns_h_on_file_lsm) )
              ALLOCATE( tmp_walltype_h_2d%var_2d(nzb_soil:nzt_soil+1,           &
-                                                1:ns_h_on_file_lsm) )          
+                                                1:ns_h_on_file_lsm) )
              ALLOCATE( tmp_walltype_h_2d2%var_2d(nzb_soil:nzt_soil,            &
-                       1:ns_h_on_file_lsm)  )                                  
-                                                                               
-          ENDIF                                                                
-                                                                               
-       CASE ( 'ns_v_on_file_lsm' )                                             
-          IF ( k == 1 )  THEN                                                  
-             READ ( 13 ) ns_v_on_file_lsm                                      
-                                                                               
-             DO  l = 0, 3                                                      
+                       1:ns_h_on_file_lsm)  )
+
+          ENDIF
+
+       CASE ( 'ns_v_on_file_lsm' )
+          IF ( k == 1 )  THEN
+             READ ( 13 ) ns_v_on_file_lsm
+
+             DO  l = 0, 3
                 IF ( ALLOCATED( tmp_walltype_v_1d(l)%var_1d ) )                &
-                   DEALLOCATE( tmp_walltype_v_1d(l)%var_1d )                   
+                   DEALLOCATE( tmp_walltype_v_1d(l)%var_1d )
                 IF ( ALLOCATED( tmp_walltype_v_2d(l)%var_2d ) )                &
-                   DEALLOCATE( tmp_walltype_v_2d(l)%var_2d )                   
+                   DEALLOCATE( tmp_walltype_v_2d(l)%var_2d )
                 IF ( ALLOCATED( tmp_walltype_v_2d2(l)%var_2d ) )               &
-                   DEALLOCATE( tmp_walltype_v_2d2(l)%var_2d )                  
-             ENDDO                                                             
-                                                                               
-!                                                                              
-!--          Allocate temporary arrays to store surface data                   
-             DO  l = 0, 3                                                      
+                   DEALLOCATE( tmp_walltype_v_2d2(l)%var_2d )
+             ENDDO
+
+!
+!--          Allocate temporary arrays to store surface data
+             DO  l = 0, 3
                 ALLOCATE( tmp_walltype_v_1d(l)                                 &
-                             %var_1d(1:ns_v_on_file_lsm(l)) )                  
+                             %var_1d(1:ns_v_on_file_lsm(l)) )
                 ALLOCATE( tmp_walltype_v_2d(l)                                 &
                              %var_2d(nzb_soil:nzt_soil+1,                      &
-                                     1:ns_v_on_file_lsm(l)) )                  
+                                     1:ns_v_on_file_lsm(l)) )
                 ALLOCATE( tmp_walltype_v_2d2(l)                                &
                              %var_2d(nzb_soil:nzt_soil,                        &
-                                     1:ns_v_on_file_lsm(l))  )                 
-             ENDDO                                                             
-                                                                               
-          ENDIF                                                                
-                                                                               
-                                                                               
-       CASE ( 'c_liq_av' )                                                     
-          IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN                             
-             ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )                         
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+                                     1:ns_v_on_file_lsm(l))  )
+             ENDDO
+
+          ENDIF
+
+
+       CASE ( 'c_liq_av' )
+          IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
+             ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                  &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-                                                                               
-       CASE ( 'c_soil_av' )                                                    
-          IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN                            
-             ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )                        
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'c_soil_av' )
+          IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
+             ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                 &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-                                                                               
-       CASE ( 'c_veg_av' )                                                     
-          IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN                             
-             ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )                         
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'c_veg_av' )
+          IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
+             ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                  &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-                                                                               
-       CASE ( 'lai_av' )                                                       
-          IF ( .NOT. ALLOCATED( lai_av ) )  THEN                               
-             ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )                           
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'lai_av' )
+          IF ( .NOT. ALLOCATED( lai_av ) )  THEN
+             ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                    &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-                                                                               
-       CASE ( 'm_liq_av' )                                                     
-          IF ( .NOT. ALLOCATED( m_liq_av ) )  THEN                             
-             ALLOCATE( m_liq_av(nysg:nyng,nxlg:nxrg) )                         
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'm_liq_av' )
+          IF ( .NOT. ALLOCATED( m_liq_av ) )  THEN
+             ALLOCATE( m_liq_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           m_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                  &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-                                                                               
-       CASE ( 'm_soil_av' )                                                    
-          IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN                            
-             ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )      
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_3d(:,:,:)                            
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'm_soil_av' )
+          IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
+             ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_3d(:,:,:)
           m_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =               &
              tmp_3d(nzb_soil:nzt_soil,nysf                                     &
-                    -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                       
-                                                                               
-       CASE ( 'qsws_liq_av' )                                                  
-          IF ( .NOT. ALLOCATED( qsws_liq_av ) )  THEN                          
-             ALLOCATE( qsws_liq_av(nysg:nyng,nxlg:nxrg) )                      
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+                    -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'qsws_liq_av' )
+          IF ( .NOT. ALLOCATED( qsws_liq_av ) )  THEN
+             ALLOCATE( qsws_liq_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           qsws_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =              &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-       CASE ( 'qsws_soil_av' )                                                 
-          IF ( .NOT. ALLOCATED( qsws_soil_av ) )  THEN                         
-             ALLOCATE( qsws_soil_av(nysg:nyng,nxlg:nxrg) )                     
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+       CASE ( 'qsws_soil_av' )
+          IF ( .NOT. ALLOCATED( qsws_soil_av ) )  THEN
+             ALLOCATE( qsws_soil_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           qsws_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =             &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-                                                                               
-       CASE ( 'qsws_veg_av' )                                                  
-          IF ( .NOT. ALLOCATED( qsws_veg_av ) )  THEN                          
-             ALLOCATE( qsws_veg_av(nysg:nyng,nxlg:nxrg) )                      
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_2d                                   
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'qsws_veg_av' )
+          IF ( .NOT. ALLOCATED( qsws_veg_av ) )  THEN
+             ALLOCATE( qsws_veg_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
           qsws_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =              &
-             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                   
-                                                                               
-       CASE ( 't_soil_av' )                                                    
-          IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN                            
-             ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )      
-          ENDIF                                                                
-          IF ( k == 1 )  READ ( 13 )  tmp_3d(:,:,:)                            
+             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 't_soil_av' )
+          IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
+             ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_3d(:,:,:)
           t_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =               &
-             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)                 
-                                                                               
-       CASE ( 'lsm_start_index_h', 'lsm_start_index_v'  )                      
-            IF ( k == 1 )  THEN                                                
-                                                                               
+             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'lsm_start_index_h', 'lsm_start_index_v'  )
+            IF ( k == 1 )  THEN
+
                IF ( ALLOCATED( start_index_on_file ) )                         &
-                  DEALLOCATE( start_index_on_file )                            
-                                                                               
+                  DEALLOCATE( start_index_on_file )
+
                ALLOCATE ( start_index_on_file(nys_on_file:nyn_on_file,         &
-               nxl_on_file:nxr_on_file) )                                      
-                                                                               
-               READ ( 13 )  start_index_on_file                                
-                                                                               
-            ENDIF                                                              
-                                                                               
-       CASE ( 'lsm_end_index_h', 'lsm_end_index_v' )                           
-            IF ( k == 1 )  THEN                                                
-                                                                               
+               nxl_on_file:nxr_on_file) )
+
+               READ ( 13 )  start_index_on_file
+
+            ENDIF
+
+       CASE ( 'lsm_end_index_h', 'lsm_end_index_v' )
+            IF ( k == 1 )  THEN
+
                IF ( ALLOCATED( end_index_on_file ) )                           &
-                  DEALLOCATE( end_index_on_file )                              
-                                                                               
+                  DEALLOCATE( end_index_on_file )
+
                ALLOCATE ( end_index_on_file(nys_on_file:nyn_on_file,           &
-                  nxl_on_file:nxr_on_file) )                                   
-                                                                               
-               READ ( 13 )  end_index_on_file                                  
-                                                                               
-            ENDIF                                                              
-                                                                               
-       CASE ( 't_soil_h' )                                                     
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                  nxl_on_file:nxr_on_file) )
+
+               READ ( 13 )  end_index_on_file
+
+            ENDIF
+
+       CASE ( 't_soil_h' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_soil_h%var_2d ) )                        &
                 ALLOCATE( t_soil_h%var_2d(nzb_soil:nzt_soil+1,                 &
-                                          1:surf_lsm_h%ns) )                   
-             READ ( 13 )  tmp_walltype_h_2d%var_2d                             
-          ENDIF                                                                
+                                          1:surf_lsm_h%ns) )
+             READ ( 13 )  tmp_walltype_h_2d%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                      t_soil_h%var_2d,                          &
                                      tmp_walltype_h_2d%var_2d,                 &
-                                     surf_lsm_h%start_index,                   & 
+                                     surf_lsm_h%start_index,                   &
                                      start_index_on_file,                      &
                                      end_index_on_file,                        &
@@ -6886 +6885 @@
                                      nxlf, nxrf, nysf, nynf,                   &
                                      nys_on_file, nyn_on_file,                 &
-                                     nxl_on_file,nxr_on_file )                 
-                                                                               
-       CASE ( 't_soil_v(0)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                     nxl_on_file,nxr_on_file )
+
+       CASE ( 't_soil_v(0)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_soil_v(0)%var_2d ) )                     &
                 ALLOCATE( t_soil_v(0)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(0)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d(0)%var_2d                          
-          ENDIF                                                                
+                                             1:surf_lsm_v(0)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d(0)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   t_soil_v(0)%var_2d,                          &
                                   tmp_walltype_v_2d(0)%var_2d,                 &
-                                  surf_lsm_v(0)%start_index,                   &  
+                                  surf_lsm_v(0)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -6905 +6904 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 't_soil_v(1)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 't_soil_v(1)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_soil_v(1)%var_2d ) )                     &
                 ALLOCATE( t_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(1)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d(1)%var_2d                          
-          ENDIF                                                                
+                                             1:surf_lsm_v(1)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d(1)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   t_soil_v(1)%var_2d,                          &
                                   tmp_walltype_v_2d(1)%var_2d,                 &
-                                  surf_lsm_v(1)%start_index,                   &   
+                                  surf_lsm_v(1)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -6924 +6923 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 't_soil_v(2)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 't_soil_v(2)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_soil_v(2)%var_2d ) )                     &
                 ALLOCATE( t_soil_v(2)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(2)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d(2)%var_2d                          
-          ENDIF                                                                
+                                             1:surf_lsm_v(2)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d(2)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   t_soil_v(2)%var_2d,                          &
                                   tmp_walltype_v_2d(2)%var_2d,                 &
-                                  surf_lsm_v(2)%start_index,                   & 
+                                  surf_lsm_v(2)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -6943 +6942 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 't_soil_v(3)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 't_soil_v(3)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_soil_v(3)%var_2d ) )                     &
                 ALLOCATE( t_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(3)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d(3)%var_2d                          
-          ENDIF                                                                
+                                             1:surf_lsm_v(3)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d(3)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   t_soil_v(3)%var_2d,                          &
                                   tmp_walltype_v_2d(3)%var_2d,                 &
-                                  surf_lsm_v(3)%start_index,                   & 
+                                  surf_lsm_v(3)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -6962 +6961 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 'm_soil_h' )                                                     
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 'm_soil_h' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_soil_h%var_2d ) )                        &
                 ALLOCATE( m_soil_h%var_2d(nzb_soil:nzt_soil+1,                 &
-                                          1:surf_lsm_h%ns) )                   
-             READ ( 13 )  tmp_walltype_h_2d2%var_2d                            
-          ENDIF                                                                
+                                          1:surf_lsm_h%ns) )
+             READ ( 13 )  tmp_walltype_h_2d2%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                     m_soil_h%var_2d,                           &
                                     tmp_walltype_h_2d2%var_2d,                 &
-                                    surf_lsm_h%start_index,                    &  
+                                    surf_lsm_h%start_index,                    &
                                     start_index_on_file,                       &
                                     end_index_on_file,                         &
@@ -6981 +6980 @@
                                     nxlf, nxrf, nysf, nynf,                    &
                                     nys_on_file, nyn_on_file,                  &
-                                    nxl_on_file,nxr_on_file )                  
-                                                                               
-       CASE ( 'm_soil_v(0)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                    nxl_on_file,nxr_on_file )
+
+       CASE ( 'm_soil_v(0)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_soil_v(0)%var_2d ) )                     &
                 ALLOCATE( m_soil_v(0)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(0)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d2(0)%var_2d                         
-          ENDIF                                                                
+                                             1:surf_lsm_v(0)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d2(0)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
-                                 m_soil_v(0)%var_2d,                           & 
+                                 m_soil_v(0)%var_2d,                           &
                                  tmp_walltype_v_2d2(0)%var_2d,                 &
-                                 surf_lsm_v(0)%start_index,                    & 
+                                 surf_lsm_v(0)%start_index,                    &
                                  start_index_on_file,                          &
                                  end_index_on_file,                            &
@@ -7000 +6999 @@
                                  nxlf, nxrf, nysf, nynf,                       &
                                  nys_on_file, nyn_on_file,                     &
-                                 nxl_on_file,nxr_on_file )                     
-                                                                               
-       CASE ( 'm_soil_v(1)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                 nxl_on_file,nxr_on_file )
+
+       CASE ( 'm_soil_v(1)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_soil_v(1)%var_2d ) )                     &
                 ALLOCATE( m_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(1)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d2(1)%var_2d                         
-          ENDIF                                                                
+                                             1:surf_lsm_v(1)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d2(1)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
-                                 m_soil_v(1)%var_2d,                           &    
+                                 m_soil_v(1)%var_2d,                           &
                                  tmp_walltype_v_2d2(1)%var_2d,                 &
-                                 surf_lsm_v(1)%start_index,                    &  
+                                 surf_lsm_v(1)%start_index,                    &
                                  start_index_on_file,                          &
                                  end_index_on_file,                            &
@@ -7019 +7018 @@
                                  nxlf, nxrf, nysf, nynf,                       &
                                  nys_on_file, nyn_on_file,                     &
-                                 nxl_on_file,nxr_on_file )                     
-                                                                               
-                                                                               
-       CASE ( 'm_soil_v(2)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                 nxl_on_file,nxr_on_file )
+
+
+       CASE ( 'm_soil_v(2)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_soil_v(2)%var_2d ) )                     &
                 ALLOCATE( m_soil_v(2)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(2)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d2(2)%var_2d                         
-          ENDIF                                                                
+                                             1:surf_lsm_v(2)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d2(2)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
-                                 m_soil_v(2)%var_2d,                           & 
+                                 m_soil_v(2)%var_2d,                           &
                                  tmp_walltype_v_2d2(2)%var_2d,                 &
-                                 surf_lsm_v(2)%start_index,                    &   
+                                 surf_lsm_v(2)%start_index,                    &
                                  start_index_on_file,                          &
                                  end_index_on_file,                            &
@@ -7039 +7038 @@
                                  nxlf, nxrf, nysf, nynf,                       &
                                  nys_on_file, nyn_on_file,                     &
-                                 nxl_on_file,nxr_on_file )                     
-                                                                               
-                                                                               
-       CASE ( 'm_soil_v(3)' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                 nxl_on_file,nxr_on_file )
+
+
+       CASE ( 'm_soil_v(3)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_soil_v(3)%var_2d ) )                     &
                 ALLOCATE( m_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,              &
-                                             1:surf_lsm_v(3)%ns) )             
-             READ ( 13 )  tmp_walltype_v_2d2(3)%var_2d                         
-          ENDIF                                                                
+                                             1:surf_lsm_v(3)%ns) )
+             READ ( 13 )  tmp_walltype_v_2d2(3)%var_2d
+          ENDIF
           CALL surface_restore_elements(                                       &
-                                 m_soil_v(3)%var_2d,                           & 
+                                 m_soil_v(3)%var_2d,                           &
                                  tmp_walltype_v_2d2(3)%var_2d,                 &
-                                 surf_lsm_v(3)%start_index,                    &  
+                                 surf_lsm_v(3)%start_index,                    &
                                  start_index_on_file,                          &
                                  end_index_on_file,                            &
@@ -7059 +7058 @@
                                  nxlf, nxrf, nysf, nynf,                       &
                                  nys_on_file, nyn_on_file,                     &
-                                 nxl_on_file,nxr_on_file )                     
-                                                                               
-                                                                               
-       CASE ( 'm_liq_h' )                                                      
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                 nxl_on_file,nxr_on_file )
+
+
+       CASE ( 'm_liq_h' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_liq_h%var_1d ) )                         &
-                ALLOCATE( m_liq_h%var_1d(1:surf_lsm_h%ns) )                    
-             READ ( 13 )  tmp_walltype_h_1d%var_1d                             
-          ENDIF                                                                
+                ALLOCATE( m_liq_h%var_1d(1:surf_lsm_h%ns) )
+             READ ( 13 )  tmp_walltype_h_1d%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                      m_liq_h%var_1d,                           &
                                      tmp_walltype_h_1d%var_1d,                 &
-                                     surf_lsm_h%start_index,                   &  
+                                     surf_lsm_h%start_index,                   &
                                      start_index_on_file,                      &
                                      end_index_on_file,                        &
@@ -7078 +7077 @@
                                      nxlf, nxrf, nysf, nynf,                   &
                                      nys_on_file, nyn_on_file,                 &
-                                     nxl_on_file,nxr_on_file )                 
-                                                                               
-                                                                               
-       CASE ( 'm_liq_v(0)' )                                                   
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                     nxl_on_file,nxr_on_file )
+
+
+       CASE ( 'm_liq_v(0)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_liq_v(0)%var_1d ) )                      &
-                ALLOCATE( m_liq_v(0)%var_1d(1:surf_lsm_v(0)%ns) )              
-             READ ( 13 )  tmp_walltype_v_1d(0)%var_1d                          
-          ENDIF                                                                
+                ALLOCATE( m_liq_v(0)%var_1d(1:surf_lsm_v(0)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(0)%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   m_liq_v(0)%var_1d,                           &
                                   tmp_walltype_v_1d(0)%var_1d,                 &
-                                  surf_lsm_v(0)%start_index,                   & 
+                                  surf_lsm_v(0)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -7097 +7096 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-                                                                               
-       CASE ( 'm_liq_v(1)' )                                                   
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+
+       CASE ( 'm_liq_v(1)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_liq_v(1)%var_1d ) )                      &
-                ALLOCATE( m_liq_v(1)%var_1d(1:surf_lsm_v(1)%ns) )              
-             READ ( 13 )  tmp_walltype_v_1d(1)%var_1d                          
-          ENDIF                                                                
+                ALLOCATE( m_liq_v(1)%var_1d(1:surf_lsm_v(1)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(1)%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   m_liq_v(1)%var_1d,                           &
                                   tmp_walltype_v_1d(1)%var_1d,                 &
-                                  surf_lsm_v(1)%start_index,                   & 
+                                  surf_lsm_v(1)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -7116 +7115 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-                                                                               
-       CASE ( 'm_liq_v(2)' )                                                   
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+
+       CASE ( 'm_liq_v(2)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_liq_v(2)%var_1d ) )                      &
-                ALLOCATE( m_liq_v(2)%var_1d(1:surf_lsm_v(2)%ns) )              
-             READ ( 13 )  tmp_walltype_v_1d(2)%var_1d                          
-          ENDIF                                                                
+                ALLOCATE( m_liq_v(2)%var_1d(1:surf_lsm_v(2)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(2)%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   m_liq_v(2)%var_1d,                           &
                                   tmp_walltype_v_1d(2)%var_1d,                 &
-                                  surf_lsm_v(2)%start_index,                   & 
+                                  surf_lsm_v(2)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -7135 +7134 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 'm_liq_v(3)' )                                                   
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 'm_liq_v(3)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( m_liq_v(3)%var_1d ) )                      &
-                ALLOCATE( m_liq_v(3)%var_1d(1:surf_lsm_v(3)%ns) )              
-             READ ( 13 )  tmp_walltype_v_1d(3)%var_1d                          
-          ENDIF                                                                
+                ALLOCATE( m_liq_v(3)%var_1d(1:surf_lsm_v(3)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(3)%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   m_liq_v(3)%var_1d,                           &
                                   tmp_walltype_v_1d(3)%var_1d,                 &
-                                  surf_lsm_v(3)%start_index,                   & 
+                                  surf_lsm_v(3)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -7153 +7152 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-                                                                               
-       CASE ( 't_surface_h' )                                                  
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+
+       CASE ( 't_surface_h' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_surface_h%var_1d ) )                     &
-                ALLOCATE( t_surface_h%var_1d(1:surf_lsm_h%ns) )                
-             READ ( 13 )  tmp_walltype_h_1d%var_1d                             
-          ENDIF                                                                
+                ALLOCATE( t_surface_h%var_1d(1:surf_lsm_h%ns) )
+             READ ( 13 )  tmp_walltype_h_1d%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                      t_surface_h%var_1d,                       &
                                      tmp_walltype_h_1d%var_1d,                 &
-                                     surf_lsm_h%start_index,                   & 
+                                     surf_lsm_h%start_index,                   &
                                      start_index_on_file,                      &
                                      end_index_on_file,                        &
@@ -7172 +7171 @@
                                      nxlf, nxrf, nysf, nynf,                   &
                                      nys_on_file, nyn_on_file,                 &
-                                     nxl_on_file,nxr_on_file )                 
-                                                                               
-       CASE ( 't_surface_v(0)' )                                               
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                     nxl_on_file,nxr_on_file )
+
+       CASE ( 't_surface_v(0)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_surface_v(0)%var_1d ) )                  &
-                ALLOCATE( t_surface_v(0)%var_1d(1:surf_lsm_v(0)%ns) )          
-             READ ( 13 )  tmp_walltype_v_1d(0)%var_1d                          
-          ENDIF                                                                
+                ALLOCATE( t_surface_v(0)%var_1d(1:surf_lsm_v(0)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(0)%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   t_surface_v(0)%var_1d,                       &
                                   tmp_walltype_v_1d(0)%var_1d,                 &
-                                  surf_lsm_v(0)%start_index,                   & 
+                                  surf_lsm_v(0)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -7190 +7189 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 't_surface_v(1)' )                                               
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 't_surface_v(1)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_surface_v(1)%var_1d ) )                  &
-                ALLOCATE( t_surface_v(1)%var_1d(1:surf_lsm_v(1)%ns) )          
-             READ ( 13 )  tmp_walltype_v_1d(1)%var_1d                          
-          ENDIF                                                                
+                ALLOCATE( t_surface_v(1)%var_1d(1:surf_lsm_v(1)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(1)%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   t_surface_v(1)%var_1d,                       &
                                   tmp_walltype_v_1d(1)%var_1d,                 &
-                                  surf_lsm_v(1)%start_index,                   & 
-                                  start_index_on_file,                         &
-                                  end_index_on_file,                           &
-                                  nxlc, nysc,                                  &
-                                  nxlf, nxrf, nysf, nynf,                      &
-                                  nys_on_file, nyn_on_file,                    & 
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 't_surface_v(2)' )                                               
-                                                                               
-          IF ( k == 1 )  THEN                                                  
-             IF ( .NOT.  ALLOCATED( t_surface_v(2)%var_1d ) )                  &
-                ALLOCATE( t_surface_v(2)%var_1d(1:surf_lsm_v(2)%ns) )          
-             READ ( 13 )  tmp_walltype_v_1d(2)%var_1d                          
-          ENDIF                                                                
-          CALL surface_restore_elements(                                       &
-                                  t_surface_v(2)%var_1d,                       &
-                                  tmp_walltype_v_1d(2)%var_1d,                 &
-                                  surf_lsm_v(2)%start_index,                   & 
+                                  surf_lsm_v(1)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -7226 +7207 @@
                                   nxlf, nxrf, nysf, nynf,                      &
                                   nys_on_file, nyn_on_file,                    &
-                                  nxl_on_file,nxr_on_file )                    
-                                                                               
-       CASE ( 't_surface_v(3)' )                                               
-                                                                               
-          IF ( k == 1 )  THEN                                                  
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 't_surface_v(2)' )
+
+          IF ( k == 1 )  THEN
+             IF ( .NOT.  ALLOCATED( t_surface_v(2)%var_1d ) )                  &
+                ALLOCATE( t_surface_v(2)%var_1d(1:surf_lsm_v(2)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(2)%var_1d
+          ENDIF
+          CALL surface_restore_elements(                                       &
+                                  t_surface_v(2)%var_1d,                       &
+                                  tmp_walltype_v_1d(2)%var_1d,                 &
+                                  surf_lsm_v(2)%start_index,                   &
+                                  start_index_on_file,                         &
+                                  end_index_on_file,                           &
+                                  nxlc, nysc,                                  &
+                                  nxlf, nxrf, nysf, nynf,                      &
+                                  nys_on_file, nyn_on_file,                    &
+                                  nxl_on_file,nxr_on_file )
+
+       CASE ( 't_surface_v(3)' )
+
+          IF ( k == 1 )  THEN
              IF ( .NOT.  ALLOCATED( t_surface_v(3)%var_1d ) )                  &
-                ALLOCATE( t_surface_v(3)%var_1d(1:surf_lsm_v(3)%ns) )          
-             READ ( 13 )  tmp_walltype_v_1d(3)%var_1d                          
-          ENDIF                                                                
+                ALLOCATE( t_surface_v(3)%var_1d(1:surf_lsm_v(3)%ns) )
+             READ ( 13 )  tmp_walltype_v_1d(3)%var_1d
+          ENDIF
           CALL surface_restore_elements(                                       &
                                   t_surface_v(3)%var_1d,                       &
                                   tmp_walltype_v_1d(3)%var_1d,                 &
-                                  surf_lsm_v(3)%start_index,                   &   
+                                  surf_lsm_v(3)%start_index,                   &
                                   start_index_on_file,                         &
                                   end_index_on_file,                           &
@@ -7272 +7271 @@
        INTEGER(iwp) ::  j       !< running index
        INTEGER(iwp) ::  m       !< running index
-       
+
        LOGICAL      ::  flag_exceed_z0  = .FALSE. !< dummy flag to indicate whether roughness length is too high
        LOGICAL      ::  flag_exceed_z0h = .FALSE. !< dummy flag to indicate whether roughness length for scalars is too high
@@ -7283 +7282 @@
        DO  m = 1, surf_lsm_h%ns
 
-          i   = surf_lsm_h%i(m)            
+          i   = surf_lsm_h%i(m)
           j   = surf_lsm_h%j(m)
-          
+
           IF ( surf_lsm_h%water_surface(m) )  THEN
 
 !
-!--          Disabled: FLake parameterization. Ideally, the Charnock 
+!--          Disabled: FLake parameterization. Ideally, the Charnock
 !--          coefficient should depend on the water depth and the fetch
 !--          length
 !             re_0 = z0(j,i) * us(j,i) / molecular_viscosity
-!        
+!
 !             z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
 !                           alpha_ch * us(j,i) / g )
@@ -7337 +7336 @@
                 flag_exceed_z0h   = .TRUE.
              ENDIF
-  
-                                 
+
+
           ENDIF
        ENDDO
@@ -7363 +7362 @@
 
     END SUBROUTINE calc_z0_water_surface
-    
+
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!>  Vertical interpolation and extrapolation of 1D soil profile input from 
-!>  dynamic input file onto the numeric vertical soil grid. 
+!>  Vertical interpolation and extrapolation of 1D soil profile input from
+!>  dynamic input file onto the numeric vertical soil grid.
 !------------------------------------------------------------------------------!
     SUBROUTINE interpolate_soil_profile( var, var_file, z_grid, z_file,        &
@@ -7452 +7451 @@
 
    END FUNCTION psi_h
-    
+
  END MODULE land_surface_model_mod