Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2608)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2696)
@@ -1,5 +1,5 @@
 !> @file land_surface_model_mod.f90
 !------------------------------------------------------------------------------!
-! This file is part of PALM.
+! This file is part of the PALM model system.
 !
 ! PALM is free software: you can redistribute it and/or modify it under the
@@ -25,4 +25,24 @@
 ! -----------------
 ! $Id$
+! Bugfix: missing USE statement for calc_mean_profile
+! do not write surface temperatures onto pt array as this might cause 
+! problems with nesting (MS)
+! Revised calculation of pt1 and qv1 (now done in surface_layer_fluxes). Bugfix
+! in calculation of surface density (cannot be done via an surface non-air 
+! temperature) (BM)
+! Bugfix: g_d was NaN for non-vegetaed surface types (BM)
+! Bugfix initialization of c_veg and lai
+! Revise data output to enable _FillValues
+! Bugfix in calcultion of r_a and rad_net_l (MS)
+! Bugfix: rad_net is not updated in case of radiation_interaction and must thu
+! be calculated again from the radiative fluxes
+! Temporary fix for cases where no soil model is used on some PEs (BM)
+! Revised input and initialization of soil and surface paramters
+! pavement_depth is variable for each surface element
+! radiation quantities belong to surface type now
+! surface fractions initialized
+! Rename lsm_last_actions into lsm_write_restart_data (MS)
+! 
+! 2608 2017-11-13 14:04:26Z schwenkel
 ! Calculation of magnus equation in external module (diagnostic_quantities_mod).
 ! Adjust calculation of vapor pressure and saturation mixing ratio that it is
@@ -265,5 +285,4 @@
 !> DALES and UCLA-LES models.
 !>
-!> @todo Restart data for vertical natural land-surfaces
 !> @todo Extensive verification energy-balance solver for vertical surfaces, 
 !>       e.g. parametrization of r_a
@@ -282,4 +301,7 @@
 !> @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 
+!>       more levels are used this may cause an potential bug
+!> @todo Routine calc_q_surface required?
 !------------------------------------------------------------------------------!
  MODULE land_surface_model_mod
@@ -287,4 +309,7 @@
     USE arrays_3d,                                                             &
         ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w
+
+    USE calc_mean_profile_mod,                                                 &
+        ONLY:  calc_mean_profile
 
     USE cloud_parameters,                                                      &
@@ -301,5 +326,13 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb 
+        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
+
+    USE netcdf_data_input_mod,                                                 &
+        ONLY :  building_type_f, init_3d, input_pids_static,                   &
+                netcdf_data_input_interpolate,                                 &
+                pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,  &
+                root_area_density_lsm_f, soil_pars_f, soil_type_f,             &
+                surface_fraction_f, vegetation_pars_f, vegetation_type_f,      &
+                water_pars_f, water_type_f
 
     USE kinds
@@ -308,7 +341,6 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  albedo, albedo_type, emissivity, force_radiation_call, rad_net, &
-               rad_sw_in, rad_lw_out, rad_lw_out_change_0, radiation_scheme,   &
-               unscheduled_radiation_calls
+        ONLY:  albedo, albedo_type, emissivity, force_radiation_call,          &
+               radiation_scheme, unscheduled_radiation_calls
         
     USE statistics,                                                            &
@@ -316,11 +348,11 @@
 
     USE surface_mod,                                                           &
-        ONLY :  surf_lsm_h, surf_lsm_v, surf_type
+        ONLY :  surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements
 
     IMPLICIT NONE
 
     TYPE surf_type_lsm
-       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1D !< 1D prognostic variable
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2D !< 2D prognostic variable
+       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1d !< 1D prognostic variable
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d !< 2D prognostic variable
     END TYPE surf_type_lsm
 
@@ -355,5 +387,5 @@
 !
 !-- LSM variables
-    CHARACTER(10) :: surface_type = 'vegetation'  !< general classification. Allowed are:
+    CHARACTER(10) :: surface_type = 'netcdf'      !< general classification. Allowed are:
                                                   !< 'vegetation', 'pavement', ('building'), 
                                                   !< 'water', and 'netcdf'
@@ -447,4 +479,5 @@
                                      m_soil_v,       & !< Soil moisture (m3/m3), vertical surface elements
                                      m_soil_v_p        !< Prog. soil moisture (m3/m3), vertical surface elements
+
 #else
     TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
@@ -474,6 +507,6 @@
     TYPE(surf_type_lsm), TARGET   ::  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 
+                                      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), DIMENSION(0:3), TARGET   ::  &
@@ -533,14 +566,14 @@
 !-- 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
-              ghf_av,       & !< average of ghf
-              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_a_av,       & !< average of r_a
-              r_s_av          !< average of r_s
+              c_liq_av,         & !< average of c_liq
+              c_soil_av,        & !< average of c_soil
+              c_veg_av,         & !< average of c_veg
+              ghf_av,           & !< average of ghf
+              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_a_av,           & !< average of r_a
+              r_s_av              !< average of r_s
                    
 
@@ -548,23 +581,23 @@
 !-- Predefined Land surface classes (vegetation_type)
     CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
-                                   'user defined              ', & !  0 
-                                   'bare soil                 ', & !  1                           
-                                   'crops, mixed farming      ', & !  2
-                                   'short grass               ', & !  3
-                                   'evergreen needleleaf trees', & !  4
-                                   'deciduous needleleaf trees', & !  5
-                                   'evergreen broadleaf trees ', & !  6
-                                   'deciduous broadleaf trees ', & !  7
-                                   'tall grass                ', & !  8
-                                   'desert                    ', & !  9
-                                   'tundra                    ', & ! 10
-                                   'irrigated crops           ', & ! 11
-                                   'semidesert                ', & ! 12
-                                   'ice caps and glaciers     ', & ! 13
-                                   'bogs and marshes          ', & ! 14
-                                   'evergreen shrubs          ', & ! 15
-                                   'deciduous shrubs          ', & ! 16
-                                   'mixed forest/woodland     ', & ! 17
-                                   'interrupted forest        '  & ! 18
+                                   'user defined              ',           & !  0 
+                                   'bare soil                 ',           & !  1                           
+                                   'crops, mixed farming      ',           & !  2
+                                   'short grass               ',           & !  3
+                                   'evergreen needleleaf trees',           & !  4
+                                   'deciduous needleleaf trees',           & !  5
+                                   'evergreen broadleaf trees ',           & !  6
+                                   'deciduous broadleaf trees ',           & !  7
+                                   'tall grass                ',           & !  8
+                                   'desert                    ',           & !  9
+                                   'tundra                    ',           & ! 10
+                                   'irrigated crops           ',           & ! 11
+                                   'semidesert                ',           & ! 12
+                                   'ice caps and glaciers     ',           & ! 13
+                                   'bogs and marshes          ',           & ! 14
+                                   'evergreen shrubs          ',           & ! 15
+                                   'deciduous shrubs          ',           & ! 16
+                                   'mixed forest/woodland     ',           & ! 17
+                                   'interrupted forest        '            & ! 18
                                                                  /)
 
@@ -572,11 +605,11 @@
 !-- Soil model classes (soil_type)
     CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
-                                   'user defined', & ! 0 
-                                   'coarse      ', & ! 1
-                                   'medium      ', & ! 2
-                                   'medium-fine ', & ! 3
-                                   'fine        ', & ! 4
-                                   'very fine   ', & ! 5
-                                   'organic     '  & ! 6
+                                   'user defined',                  & ! 0 
+                                   'coarse      ',                  & ! 1
+                                   'medium      ',                  & ! 2
+                                   'medium-fine ',                  & ! 3
+                                   'fine        ',                  & ! 4
+                                   'very fine   ',                  & ! 5
+                                   'organic     '                   & ! 6
                                                                  /)
 
@@ -600,24 +633,51 @@
                                    'artifical turf (sports)      ', & ! 14
                                    'clay (sports)                '  & ! 15
-                                                                 /)                                                                 
-   
-
-
-
-                                                              
+                                                                 /)                                                              
+                                                                 
 !
 !-- Water classes
     CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
-                                   'user defined', & ! 0 
-                                   'lake        ', & ! 1
-                                   'river       ', & ! 2
-                                   'ocean       ', & ! 3
-                                   'pond        ', & ! 4
-                                   'fountain    '  & ! 5
+                                   'user defined',                   & ! 0 
+                                   'lake        ',                   & ! 1
+                                   'river       ',                   & ! 2
+                                   'ocean       ',                   & ! 3
+                                   'pond        ',                   & ! 4
+                                   'fountain    '                    & ! 5
                                                                   /)                                                                  
-                                                                 
-                                                                 !
+!
+!-- IDs for vegetation, pavement and water surfaces
+    INTEGER(iwp) ::  ind_veg = 0    !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc.      
+    INTEGER(iwp) ::  ind_pav = 1    !< index for pavement surfaces, used to assess surface-fraction, albedo, etc.        
+    INTEGER(iwp) ::  ind_wat = 2    !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc.                          
+                    
+!
 !-- Land surface parameters according to the respective classes (vegetation_type)
-
+    INTEGER(iwp) ::  ind_v_rc_min = 0    !< index for r_canopy_min in vegetation_pars
+    INTEGER(iwp) ::  ind_v_rc_lai = 1    !< index for LAI in vegetation_pars
+    INTEGER(iwp) ::  ind_v_c_veg   = 2   !< index for c_veg in vegetation_pars
+    INTEGER(iwp) ::  ind_v_gd  = 3       !< index for g_d in vegetation_pars
+    INTEGER(iwp) ::  ind_v_z0 = 4        !< index for z0 in vegetation_pars
+    INTEGER(iwp) ::  ind_v_z0qh = 5      !< index for z0h / z0q in vegetation_pars
+    INTEGER(iwp) ::  ind_v_lambda_s = 6  !< index for lambda_s_s in vegetation_pars
+    INTEGER(iwp) ::  ind_v_lambda_u = 7  !< index for lambda_s_u in vegetation_pars
+    INTEGER(iwp) ::  ind_v_f_sw_in = 8   !< index for f_sw_in in vegetation_pars
+    INTEGER(iwp) ::  ind_v_c_surf = 9    !< index for c_surface in vegetation_pars
+    INTEGER(iwp) ::  ind_v_at = 10       !< index for albedo_type in vegetation_pars
+    INTEGER(iwp) ::  ind_v_emis = 11     !< index for emissivity in vegetation_pars
+
+    INTEGER(iwp) ::  ind_w_temp     = 0    !< index for temperature in water_pars
+    INTEGER(iwp) ::  ind_w_z0       = 1    !< index for z0 in water_pars
+    INTEGER(iwp) ::  ind_w_z0h      = 2    !< index for z0h in water_pars
+    INTEGER(iwp) ::  ind_w_lambda_s = 3    !< index for lambda_s_s in water_pars
+    INTEGER(iwp) ::  ind_w_lambda_u = 4    !< index for lambda_s_u in water_pars
+    INTEGER(iwp) ::  ind_w_at       = 5    !< index for albedo type in water_pars
+    INTEGER(iwp) ::  ind_w_emis     = 6    !< index for emissivity in water_pars
+
+    INTEGER(iwp) ::  ind_p_z0       = 0    !< index for z0 in pavement_pars
+    INTEGER(iwp) ::  ind_p_z0h      = 1    !< index for z0h in pavement_pars
+    INTEGER(iwp) ::  ind_p_at       = 2    !< index for albedo type in pavement_pars
+    INTEGER(iwp) ::  ind_p_emis     = 3    !< index for emissivity in pavement_pars
+    INTEGER(iwp) ::  ind_p_lambda_h = 0    !< index for lambda_h in pavement_subsurface_pars
+    INTEGER(iwp) ::  ind_p_rho_c    = 1    !< index for rho_c in pavement_pars
 !
 !-- Land surface parameters
@@ -649,22 +709,22 @@
 !--                                level 1 - level 4 according to zs_ref
     REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
-                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & !  1
-                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & !  2
-                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & !  3
-                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & !  4
-                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & !  5
-                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & !  6
-                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & !  7
-                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & !  8
-                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & !  9
-                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 10
-                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 11
-                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 12
-                                 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.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
+                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  1
+                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  2
+                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  3
+                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  4
+                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  5
+                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  6
+                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  7
+                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  8
+                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  9
+                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & ! 10
+                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 11
+                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 12
+                                 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.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
                                  /), (/ 4, 18 /) )
 
@@ -682,4 +742,5 @@
                       1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6
                                                                      /), (/ 8, 6 /) )
+
 
 !
@@ -747,12 +808,12 @@
 !
 !-- TO BE FILLED
-!-- Water parameters  temperature,     z0,      z0h, albedo_type, emissivity,
-    REAL(wp), DIMENSION(0:4,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
-                         283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 1
-                         283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 2
-                         283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 3
-                         283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 4
-                         283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp  & ! 5
-                                                                /), (/ 5, 5 /) )                                                                     
+!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
+    REAL(wp), DIMENSION(0:6,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
+       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 1
+       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 2
+       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 3
+       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 /) )                                                                   
                                                                                                                                       
     SAVE
@@ -764,9 +825,10 @@
 !
 !-- Public functions
-    PUBLIC lsm_check_data_output, lsm_check_data_output_pr,                    &
+    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_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
            lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
-           lsm_swap_timelevel, lsm_read_restart_data, lsm_last_actions
+           lsm_swap_timelevel, lsm_read_restart_data, lsm_write_restart_data
 ! !vegetat
 !-- Public parameters, constants and initial values
@@ -781,4 +843,7 @@
     PUBLIC m_soil_h, t_soil_h
 
+    INTERFACE lsm_boundary_condition
+       MODULE PROCEDURE lsm_boundary_condition
+    END INTERFACE lsm_boundary_condition
 
     INTERFACE lsm_check_data_output
@@ -842,9 +907,74 @@
     END INTERFACE lsm_read_restart_data
 
-    INTERFACE lsm_last_actions
-       MODULE PROCEDURE lsm_last_actions
-    END INTERFACE lsm_last_actions
+    INTERFACE lsm_write_restart_data
+       MODULE PROCEDURE lsm_write_restart_data
+    END INTERFACE lsm_write_restart_data
 
  CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set internal Neumann boundary condition at outer soil grid points 
+!> for temperature and humidity. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE lsm_boundary_condition
+ 
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: i      !< grid index x-direction
+    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
+    INTEGER(iwp) :: j      !< grid index y-direction
+    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
+    INTEGER(iwp) :: k      !< grid index z-direction
+    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
+    INTEGER(iwp) :: l      !< running index surface-orientation
+    INTEGER(iwp) :: m      !< running index surface elements
+
+    koff = surf_lsm_h%koff
+    DO  m = 1, surf_lsm_h%ns
+       i = surf_lsm_h%i(m)
+       j = surf_lsm_h%j(m)
+       k = surf_lsm_h%k(m)
+       pt(k+koff,j,i) = pt(k,j,i)
+    ENDDO
+
+    DO  l = 0, 3
+       ioff = surf_lsm_v(l)%ioff
+       joff = surf_lsm_v(l)%joff
+       DO  m = 1, surf_lsm_v(l)%ns
+          i = surf_lsm_v(l)%i(m)
+          j = surf_lsm_v(l)%j(m)
+          k = surf_lsm_v(l)%k(m)
+          pt(k,j+joff,i+ioff) = pt(k,j,i)
+       ENDDO
+    ENDDO
+!
+!-- In case of humidity, set boundary conditions also for q and vpt.
+    IF ( humidity )  THEN
+       koff = surf_lsm_h%koff
+       DO  m = 1, surf_lsm_h%ns
+          i = surf_lsm_h%i(m)
+          j = surf_lsm_h%j(m)
+          k = surf_lsm_h%k(m)
+          q(k+koff,j,i)   = q(k,j,i)
+          vpt(k+koff,j,i) = vpt(k,j,i)
+       ENDDO
+
+       DO  l = 0, 3
+          ioff = surf_lsm_v(l)%ioff
+          joff = surf_lsm_v(l)%joff
+          DO  m = 1, surf_lsm_v(l)%ns
+             i = surf_lsm_v(l)%i(m)
+             j = surf_lsm_v(l)%j(m)
+             k = surf_lsm_v(l)%k(m)
+             q(k,j+joff,i+ioff)   = q(k,j,i)
+             vpt(k,j+joff,i+ioff) = vpt(k,j,i)
+          ENDDO
+       ENDDO
+    ENDIF
+
+ END SUBROUTINE lsm_boundary_condition
 
 !------------------------------------------------------------------------------!
@@ -975,4 +1105,5 @@
 
  END SUBROUTINE lsm_check_data_output
+
 
 
@@ -1060,8 +1191,5 @@
         ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
                most_method
-                 
-    USE radiation_model_mod,                                                   &
-        ONLY:  radiation
-    
+                     
     
     IMPLICIT NONE
@@ -1075,8 +1203,8 @@
          TRIM( surface_type ) /= 'pavement'    .AND.                           &
          TRIM( surface_type ) /= 'water'       .AND.                           &
-         TRIM( surface_type ) /= 'netcdf' )  THEN
-          message_string = 'unknown surface type surface_type = "' //          &
+         TRIM( surface_type ) /= 'netcdf' )  THEN  
+       message_string = 'unknown surface type surface_type = "' //             &
                         TRIM( surface_type ) // '"'
-                CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
+       CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -1164,5 +1292,5 @@
           ENDIF
        ENDIF
-     
+
        IF ( vegetation_type == 1 )  THEN
           IF ( vegetation_coverage /= 9999999.9_wp  .AND.  vegetation_coverage &
@@ -1260,11 +1388,30 @@
     
     ENDIF
+
+    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
+       IF ( ANY( water_type_f%var /= water_type_f%fill )  .AND.                &
+            TRIM( most_method ) == 'lookup' )  THEN    
+          WRITE( message_string, * ) 'water-surfaces are not allowed in ' //   &
+                                     'combination with most_method = ',        &
+                                     TRIM( most_method )
+          CALL message( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
+       ENDIF
+!
+!--    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( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
+       ENDIF
+    ENDIF
     
 !
 !-- Temporary message as long as NetCDF input is not available
-    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
-             message_string = 'surface_type = netcdf '//                       &
-                              'is not supported at the moment'
-             CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
+    IF ( TRIM( surface_type ) == 'netcdf'  .AND.  .NOT.  input_pids_static )   &
+    THEN
+       message_string = 'surface_type = netcdf requires static input file.'
+       CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -1362,5 +1509,5 @@
     DO  k = nzb_soil, nzt_soil
        IF ( soil_moisture(k) > saturation_moisture )  THEN 
-          message_string = 'soil_moisture must not exceed its saturation' //   &
+          message_string = 'soil_moisture must not exceed its saturation' //    &
                             ' value'
           CALL message( 'check_parameters', 'PA0458', 1, 2, 0, 6, 0 )
@@ -1378,4 +1525,5 @@
     ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
     ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
+
 
     zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
@@ -1428,5 +1576,4 @@
     INTEGER(iwp) ::  k         !< running index
     INTEGER(iwp) ::  k_off     !< offset to determine index of surface element, seen from atmospheric grid point, for z
-    INTEGER(iwp) ::  k_rad     !< index to access radiation array
     INTEGER(iwp) ::  ks        !< running index
     INTEGER(iwp) ::  l         !< surface-facing index
@@ -1454,7 +1601,5 @@
                 lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
                 lambda_surface, & !< Current value of lambda_surface (W/m2/K)
-                m_liq_max,   & !< maxmimum value of the liq. water reservoir
-                pt1,         & !< potential temperature at first grid level
-                qv1            !< specific humidity at first grid level
+                m_liq_max      !< maxmimum value of the liq. water reservoir
 
     TYPE(surf_type_lsm), POINTER ::  surf_t_surface
@@ -1472,4 +1617,5 @@
     IF ( horizontal )  THEN
        surf              => surf_lsm_h
+
        surf_t_surface    => t_surface_h
        surf_t_surface_p  => t_surface_h_p
@@ -1480,10 +1626,7 @@
        surf_m_soil       => m_soil_h
        surf_t_soil       => t_soil_h
-
-       k_off     = -1
-       j_off     = 0
-       i_off     = 0
     ELSE
        surf              => surf_lsm_v(l)
+
        surf_t_surface    => t_surface_v(l)
        surf_t_surface_p  => t_surface_v_p(l)
@@ -1494,20 +1637,12 @@
        surf_m_soil       => m_soil_v(l)
        surf_t_soil       => t_soil_v(l)
-
-       k_off = 0
-       IF ( l == 0 )  THEN 
-          j_off = -1 
-          i_off = 0
-       ELSEIF ( l == 1 )  THEN 
-          j_off = 1 
-          i_off = 0
-       ELSEIF ( l == 2 )  THEN 
-          j_off = 0 
-          i_off = -1
-       ELSEIF ( l == 3 )  THEN 
-          j_off = 0 
-          i_off = 1
-       ENDIF
     ENDIF
+
+!
+!-- Index offset of surface element point with respect to adjoining 
+!-- atmospheric grid point 
+    k_off = surf%koff
+    j_off = surf%joff
+    i_off = surf%ioff
 
 !
@@ -1520,10 +1655,4 @@
        j   = surf%j(m)
        k   = surf%k(m)
-!
-!--    Determine height index for radiation. Note, in clear-sky case radiation
-!--    arrays have rank 0 in first dimensions, so index must be zero. In case
-!--    of RRTMG radiation arrays have non-zero rank in first dimension, so that 
-!--    radiation can be obtained at respective height level
-       k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )
 
 !
@@ -1589,18 +1718,18 @@
 !--    time step is used here. Note that this formulation is the
 !--    equivalent to the ECMWF formulation using drag coefficients
-       IF ( cloud_physics )  THEN
-          pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
-          qv1 = q(k,j,i) - ql(k,j,i)
-       ELSEIF ( cloud_droplets ) THEN
-          pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
-          qv1 = q(k,j,i) 
-       ELSE
-          pt1 = pt(k,j,i)
-          IF ( humidity )  THEN
-             qv1 = q(k,j,i)
-          ELSE
-             qv1 = 0.0_wp
-          ENDIF
-       ENDIF
+!        IF ( cloud_physics )  THEN
+!           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+!           qv1 = q(k,j,i) - ql(k,j,i)
+!        ELSEIF ( cloud_droplets ) THEN
+!           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+!           qv1 = q(k,j,i) 
+!        ELSE
+!           pt1 = pt(k,j,i)
+!           IF ( humidity )  THEN
+!              qv1 = q(k,j,i)
+!           ELSE
+!              qv1 = 0.0_wp
+!           ENDIF
+!        ENDIF
 !
 !--    Calculate aerodynamical resistance. For horizontal and vertical 
@@ -1608,19 +1737,23 @@
 !--    can be obtain via parameterization of Mason (2000) / 
 !--    Krayenhoff and Voogt (2006).
-!--    To do: detailed investigation which approach is better!
+!--    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. 
        IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
-          surf%r_a(m) = ( pt1 - surf_lsm_h%pt_surface(m) ) /                   &
-                        ( surf%ts(m) * surf%us(m) + 1.0E-20_wp )
+          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
+                             ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
        ELSE
           surf%r_a(m) = 1.0_wp / ( 11.8_wp + 4.2_wp *                          &
-                        SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +      &
-                              ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +      &
-                              ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 )      &
+                        SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
+                                   ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
+                                   ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
+                              0.01_wp ) )                                      &
                                  )
        ENDIF
 !
 !--    Make sure that the resistance does not drop to zero for neutral 
-!--    stratification
-       IF ( ABS( surf%r_a(m) ) < 1.0_wp )  surf%r_a(m) = 1.0_wp
+!--    stratification. 
+       IF ( surf%r_a(m) < 1.0_wp )  surf%r_a(m) = 1.0_wp
 !
 !--    Second step: calculate canopy resistance r_canopy
@@ -1629,7 +1762,7 @@
 !--    f1: correction for incoming shortwave radiation (stomata close at 
 !--    night)
-       f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off)         &
-                 + 0.05_wp ) / (0.81_wp * (0.004_wp                            &
-                 * rad_sw_in(k_rad,j+j_off,i+i_off) + 1.0_wp)) )
+       f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /          &
+                        (0.81_wp * (0.004_wp * surf%rad_sw_in(m)               &
+                         + 1.0_wp)) )
 !
 !--    f2: correction for soil moisture availability to plants (the 
@@ -1665,5 +1798,5 @@
 !
 !--       Calculate vapour pressure
-          e  = qv1 * surface_pressure / ( qv1 + 0.622_wp )
+          e  = surf%qv1(m) * surface_pressure / ( surf%qv1(m) + 0.622_wp )
           f3 = EXP ( - surf%g_d(m) * (e_s - e) )
        ELSE
@@ -1682,6 +1815,6 @@
 
 
-       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min )                         &
-            / ( surf%m_fc(nzb_soil,m) - m_min )
+       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) /                       &
+            ( surf%m_fc(nzb_soil,m) - m_min )
        f2 = MAX( f2, 1.0E-20_wp )
        f2 = MIN( f2, 1.0_wp     )
@@ -1698,5 +1831,5 @@
        IF ( surf%pavement_surface(m) )  THEN
           m_liq_max = m_max_depth * 5.0_wp
-          surf%c_liq(m) = MIN( 1.0_wp, (surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
+          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
        ELSE
           m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)              &
@@ -1710,5 +1843,5 @@
 !--    In case of dewfall, set evapotranspiration to zero
 !--    All super-saturated water is then removed from the air
-       IF ( humidity  .AND.  q_s <= qv1 )  THEN
+       IF ( humidity  .AND.  q_s <= surf%qv1(m) )  THEN
           surf%r_canopy(m) = 0.0_wp
           surf%r_soil(m)   = 0.0_wp
@@ -1747,6 +1880,9 @@
        dq_s_dt = 0.622_wp * e_s_dt / ( surface_pressure - e_s_dt )
 !
-!--    Add LW up so that it can be removed in prognostic equation
-       surf%rad_net_l(m) = rad_net(j,i) + rad_lw_out(nzb,j,i)
+!--    Calculate net radiation radiation without longwave outgoing flux because
+!--    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)              &
+                           + surf%rad_lw_in(m)
 !
 !--    Calculate new skin temperature
@@ -1754,25 +1890,27 @@
 !
 !--       Numerator of the prognostic equation
-          coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
-                   * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
-                   + f_shf * pt1 + f_qsws * ( qv1 - q_s                        &
+          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
+                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
+                   + f_shf * surf%pt1(m) + f_qsws * ( surf%qv1(m) - q_s        &
                    + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface     &
                    * surf_t_soil%var_2d(nzb_soil,m)
+
 !
 !--       Denominator of the prognostic equation
-          coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt                 &
+          coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt              &
                    + lambda_surface + f_shf / exn
        ELSE
 !
 !--       Numerator of the prognostic equation
-          coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
-                   * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
-                   + f_shf * pt1  + lambda_surface                             &
+          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
+                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
+                   + f_shf * surf%pt1(m)  + lambda_surface                     &
                    * surf_t_soil%var_2d(nzb_soil,m)
 !
 !--       Denominator of the prognostic equation
-          coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
+          coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exn
 
        ENDIF
+
 
        tend = 0.0_wp
@@ -1819,9 +1957,21 @@
 !--    computing time for the radiation code).
        IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )       &
-            > 0.2_wp  .AND. unscheduled_radiation_calls )  THEN
+            > 0.2_wp  .AND. &
+            unscheduled_radiation_calls )  THEN
           force_radiation_call_l = .TRUE.
        ENDIF
 
-       pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
+
+!        IF ( surf_t_surface_p%var_1d(m) < 270.0_wp .OR. surf_t_surface_p%var_1d(m) > 350.0_wp )  THEN
+!           PRINT*, "myid: ", myid
+!           PRINT*, "m: ", m
+!           PRINT*, "i,j,k: ", i, j, k
+!           PRINT*, "pt_p: ", surf_t_surface_p%var_1d(m)
+!           PRINT*, "f_shf: ", f_shf
+!           PRINT*, "f_qsws: ", f_qsws
+!         ENDIF
+
+
+!        pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
        surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exn
 
@@ -1829,10 +1979,10 @@
 !--    Calculate fluxes
        surf%rad_net_l(m) = surf%rad_net_l(m) +                                 &
-                            rad_lw_out_change_0(j,i)                           &
-                          * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)     &
-                          - rad_lw_out_change_0(j,i) * surf_t_surface_p%var_1d(m)
-
-       rad_net(j,i) = surf%rad_net_l(m)
-       rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) *  &
+                            surf%rad_lw_out_change_0(m)                        &
+                          * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)      &
+                          - surf%rad_lw_out_change_0(m) * surf_t_surface_p%var_1d(m)
+
+       surf%rad_net(m) = surf%rad_net_l(m)
+       surf%rad_lw_out(m) = surf%rad_lw_out(m) + surf%rad_lw_out_change_0(m) * &
                      ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
 
@@ -1840,21 +1990,20 @@
                                              - surf_t_soil%var_2d(nzb_soil,m) )
 
-       surf%shf(m) = - f_shf * ( pt1 - surf%pt_surface(m) ) / cp
-
+       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / cp
 
        IF ( humidity )  THEN
-          surf%qsws(m)  = - f_qsws * ( qv1 - q_s + dq_s_dt                     &
+          surf%qsws(m)  = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt             &
                           * surf_t_surface%var_1d(m) - dq_s_dt *               &
                             surf_t_surface_p%var_1d(m) )
 
-          surf%qsws_veg(m)  = - f_qsws_veg  * ( qv1 - q_s                      &
+          surf%qsws_veg(m)  = - f_qsws_veg  * ( surf%qv1(m) - q_s              &
                               + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
                               * surf_t_surface_p%var_1d(m) )
 
-          surf%qsws_soil(m) = - f_qsws_soil * ( qv1 - q_s                      &
+          surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s              &
                               + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
                               * surf_t_surface_p%var_1d(m) )
 
-          surf%qsws_liq(m)  = - f_qsws_liq  * ( qv1 - q_s                      &
+          surf%qsws_liq(m)  = - f_qsws_liq  * ( surf%qv1(m) - q_s              &
                               + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
                               * surf_t_surface_p%var_1d(m) )
@@ -1866,5 +2015,5 @@
           surf%r_s(m) = 1.0E10_wp
        ELSE
-          surf%r_s(m) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
+          surf%r_s(m) = - rho_lv * ( surf%qv1(m) - q_s + dq_s_dt               &
                           * surf_t_surface%var_1d(m) - dq_s_dt *               &
                             surf_t_surface_p%var_1d(m) ) /                     &
@@ -1954,5 +2103,5 @@
                       intermediate_timestep_count_max )  THEN
                 surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +                  &
-                                              5.3125_wp * surf_tm_liq_m%var_1d(m)
+                                           5.3125_wp * surf_tm_liq_m%var_1d(m)
              ENDIF
           ENDIF
@@ -2015,5 +2164,5 @@
           q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
 
-          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )
+          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp )
 
 !
@@ -2028,5 +2177,4 @@
                                           q(k,j,i)
           ENDIF
-
 !
 !--       Update virtual potential temperature
@@ -2098,7 +2246,11 @@
        ENDIF
 
-       WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),           &
-                        TRIM (soil_type_name(soil_type) )
-       WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
+       IF ( vegetation_type_f%from_file )  THEN
+          WRITE( io, 5 )
+       ELSE
+          WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),        &
+                         TRIM (soil_type_name(soil_type) )
+       ENDIF
+       WRITE( io, 6 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
                         TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
                         TRIM( vertical_index_chr )
@@ -2111,5 +2263,7 @@
 4   FORMAT ('    --> Land surface type  : ',A,/                                &
             '    --> Soil porosity type : ',A)
-5   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
+5   FORMAT ('    --> Land surface type  : read from file' /                    &
+            '    --> Soil porosity type : read from file' )
+6   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
             '       Height:        ',A,'  m'/                                  &
             '       Temperature:   ',A,'  K'/                                  &
@@ -2117,4 +2271,5 @@
             '       Root fraction: ',A,'  '/                                   &
             '       Grid point:    ',A)
+
 
     END SUBROUTINE lsm_header
@@ -2133,19 +2288,18 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i     !< running index
-       INTEGER(iwp) ::  i_off !< index offset of surface element, seen from atmospheric grid point 
-       INTEGER(iwp) ::  j     !< running index
-       INTEGER(iwp) ::  j_off !< index offset of surface element, seen from atmospheric grid point 
-       INTEGER(iwp) ::  k     !< running index
-       INTEGER(iwp) ::  kn    !< running index
-       INTEGER(iwp) ::  ko    !< running index
-       INTEGER(iwp) ::  kroot !< running index
-       INTEGER(iwp) ::  kzs   !< running index
-       INTEGER(iwp) ::  l     !< running index surface facing
-       INTEGER(iwp) ::  m     !< running index
+       INTEGER(iwp) ::  i                       !< running index
+       INTEGER(iwp) ::  i_off                   !< index offset of surface element, seen from atmospheric grid point 
+       INTEGER(iwp) ::  j                       !< running index
+       INTEGER(iwp) ::  j_off                   !< index offset of surface element, seen from atmospheric grid point 
+       INTEGER(iwp) ::  k                       !< running index
+       INTEGER(iwp) ::  kn                      !< running index
+       INTEGER(iwp) ::  ko                      !< running index
+       INTEGER(iwp) ::  kroot                   !< running index
+       INTEGER(iwp) ::  kzs                     !< running index
+       INTEGER(iwp) ::  l                       !< running index surface facing
+       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 
-
-
-       REAL(wp) :: pt1   !< potential temperature at first grid level
+       INTEGER(iwp) ::  n_surf                  !< number of surface types of given surface element 
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
@@ -2157,5 +2311,6 @@
 !--    If no cloud physics is used, rho_surface has not been calculated before
        IF (  .NOT.  cloud_physics )  THEN
-          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
+          CALL calc_mean_profile( pt, 4 )
+          rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(nzb+1,1,4,0) * exn )
        ENDIF
 
@@ -2206,5 +2361,18 @@
           surf_lsm_v(l)%r_soil     = 0.0_wp
        ENDDO
-   
+
+!
+!--    Set initial values for prognostic soil quantities
+       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+          t_soil_h%var_2d = 0.0_wp
+          m_soil_h%var_2d = 0.0_wp
+          m_liq_h%var_1d  = 0.0_wp
+
+          DO  l = 0, 3
+             t_soil_v(l)%var_2d = 0.0_wp
+             m_soil_v(l)%var_2d = 0.0_wp
+             m_liq_v(l)%var_1d  = 0.0_wp
+          ENDDO
+       ENDIF
 !
 !--    Allocate 3D soil model arrays
@@ -2261,50 +2429,143 @@
           ENDIF      
        ENDDO
-
+!
+!--    Allocate albedo type and emissivity for vegetation, water and pavement 
+!--    fraction.
+!--    Set default values at each surface element. 
+       ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) )
+       ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) )
+       surf_lsm_h%albedo_type = 0      
+       surf_lsm_h%emissivity  = emissivity
+       DO  l = 0, 3
+          ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) )
+          ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns)  )
+          surf_lsm_v(l)%albedo_type = 0
+          surf_lsm_v(l)%emissivity  = emissivity
+       ENDDO
+!
+!--    Allocate arrays for relative surface fraction. 
+!--    0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
+       ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) )
+       surf_lsm_h%frac = 0.0_wp
+       DO  l = 0, 3
+          ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) )
+          surf_lsm_v(l)%frac = 0.0_wp
+       ENDDO
+!
+!--    For vertical walls only - allocate special flag indicating if any building is on
+!--    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) )
+       ENDDO
 !
 !--    Set flag parameter for the prescribed surface type depending on user 
-!--    input.
-       DO  m = 1, surf_lsm_h%ns
-
-          SELECT CASE ( TRIM( surface_type ) )
+!--    input. Set surface fraction to 1 for the respective type.
+       SELECT CASE ( TRIM( surface_type ) )
           
-             CASE ( 'vegetation' )
+          CASE ( 'vegetation' )
          
-                surf_lsm_h%vegetation_surface(m)    = .TRUE.
+             surf_lsm_h%vegetation_surface = .TRUE.
+             surf_lsm_h%frac(ind_veg,:) = 1.0_wp
+             DO  l = 0, 3
+                surf_lsm_v(l)%vegetation_surface = .TRUE.
+                surf_lsm_v(l)%frac(ind_veg,:) = 1.0_wp
+             ENDDO
    
-             CASE ( 'water' )
+          CASE ( 'water' )
              
-                surf_lsm_h%water_surface(m)      = .TRUE.
-
-             CASE ( 'pavement' )
+             surf_lsm_h%water_surface = .TRUE.
+             surf_lsm_h%frac(ind_wat,:) = 1.0_wp
+!
+!--          Note, vertical water surface does not really make sense.
+             DO  l = 0, 3  
+                surf_lsm_v(l)%water_surface   = .TRUE.
+                surf_lsm_v(l)%frac(ind_wat,:) = 1.0_wp
+             ENDDO
+
+          CASE ( 'pavement' )
              
-                surf_lsm_h%pavement_surface(m)   = .TRUE.
-    
-          END SELECT
-  
-       ENDDO
-       
-       DO  l = 0, 3
-          SELECT CASE ( TRIM( surface_type ) )
-          
-             CASE ( 'vegetation' )
-         
-                surf_lsm_v(l)%vegetation_surface = .TRUE.
-   
-             CASE ( 'water' )
-
-                surf_lsm_v(l)%water_surface   = .TRUE.
-
-             CASE ( 'pavement' )
-             
-                surf_lsm_h%pavement_surface   = .TRUE.
-
-          END SELECT   
-
-       ENDDO
-!
-!--    Initialize standard soil types. It is possible to overwrite each
-!--    parameter by setting the respecticy NAMELIST variable to a 
-!--    value /= 9999999.9. 
+             surf_lsm_h%pavement_surface = .TRUE.
+                surf_lsm_h%frac(ind_pav,:) = 1.0_wp
+             DO  l = 0, 3
+                surf_lsm_v(l)%pavement_surface   = .TRUE.
+                surf_lsm_v(l)%frac(ind_pav,:) = 1.0_wp
+             ENDDO
+
+          CASE ( 'netcdf' )
+
+             DO  m = 1, surf_lsm_h%ns
+                i = surf_lsm_h%i(m)
+                j = surf_lsm_h%j(m)
+                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )        &
+                   surf_lsm_h%vegetation_surface(m) = .TRUE.
+                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )          &
+                   surf_lsm_h%pavement_surface(m) = .TRUE.
+                IF ( water_type_f%var(j,i)      /= water_type_f%fill )             &
+                   surf_lsm_h%water_surface(m) = .TRUE.
+             ENDDO
+             DO  l = 0, 3
+                DO  m = 1, surf_lsm_v(l)%ns
+!
+!--                Only for vertical surfaces. Check if natural walls at reference
+!--                grid point are covered by any building. This case, problems 
+!--                with initialization will aris if index offsets are used. 
+!--                In order to deal with this, set special flag. 
+                   surf_lsm_v(l)%building_covered(m) = .FALSE.
+                   IF ( building_type_f%from_file )  THEN
+                      i = surf_lsm_v(l)%i(m) + surf_lsm_v(l)%ioff
+                      j = surf_lsm_v(l)%j(m) + surf_lsm_v(l)%joff
+                      IF ( building_type_f%var(j,i) /= 0 )                     &
+                         surf_lsm_v(l)%building_covered(m) = .TRUE.
+                   ENDIF
+!
+!--                Normally proceed with setting surface types.
+                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,     &
+                                            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)%vegetation_surface(m) = .TRUE.
+                   IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )   &
+                      surf_lsm_v(l)%pavement_surface(m) = .TRUE.
+                   IF ( water_type_f%var(j,i)      /= water_type_f%fill )      &
+                      surf_lsm_v(l)%water_surface(m) = .TRUE.
+                ENDDO
+             ENDDO
+
+       END SELECT
+!
+!--    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 ( input_pids_static )  THEN
+          DO  m = 1, surf_lsm_h%ns
+             i = surf_lsm_h%i(m)
+             j = surf_lsm_h%j(m)
+!
+!--          0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction             
+             surf_lsm_h%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i)         
+             surf_lsm_h%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i)        
+             surf_lsm_h%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i)
+
+          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) ) 
+                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)%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i)         
+                surf_lsm_v(l)%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i)        
+                surf_lsm_v(l)%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i)
+
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    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  
  
@@ -2339,131 +2600,4 @@
           IF ( residual_moisture == 9999999.9_wp )  THEN
              residual_moisture = soil_pars(7,soil_type)       
-          ENDIF
-
-       ENDIF    
-  
-!
-!--    Check whether parameters from the lookup tables are to be used
-       IF ( vegetation_type /= 0 )  THEN
-
-          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
-             min_canopy_resistance = vegetation_pars(0,vegetation_type)
-          ENDIF
-
-          IF ( leaf_area_index == 9999999.9_wp )  THEN
-             leaf_area_index = vegetation_pars(1,vegetation_type)         
-          ENDIF
-
-          IF ( vegetation_coverage == 9999999.9_wp )  THEN
-             vegetation_coverage = vegetation_pars(2,vegetation_type)     
-          ENDIF
-
-          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
-              canopy_resistance_coefficient= vegetation_pars(3,vegetation_type)     
-          ENDIF
-
-          IF ( z0_vegetation == 9999999.9_wp )  THEN
-             z0_vegetation  = vegetation_pars(4,vegetation_type) 
-          ENDIF
-
-          IF ( z0h_vegetation == 9999999.9_wp )  THEN
-             z0h_vegetation = vegetation_pars(5,vegetation_type)
-          ENDIF   
-          
-          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-             lambda_surface_stable = vegetation_pars(6,vegetation_type) 
-          ENDIF
-
-          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-             lambda_surface_unstable = vegetation_pars(7,vegetation_type)           
-          ENDIF
-
-          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
-             f_shortwave_incoming = vegetation_pars(8,vegetation_type)        
-          ENDIF
-
-          IF ( c_surface == 9999999.9_wp )  THEN
-             c_surface = vegetation_pars(9,vegetation_type)        
-          ENDIF
- 
-          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
-             albedo_type = INT(vegetation_pars(10,vegetation_type))       
-          ENDIF
-   
-          IF ( emissivity == 9999999.9_wp )  THEN
-             emissivity = vegetation_pars(11,vegetation_type)     
-          ENDIF
-          
-       ENDIF
-
-       IF ( water_type /= 0 )  THEN
-
-          IF ( water_temperature == 9999999.9_wp )  THEN
-             water_temperature = water_pars(0,water_type)       
-          ENDIF
-
-          IF ( z0_water == 9999999.9_wp )  THEN
-             z0_water = water_pars(1,water_type)
-          ENDIF       
-
-          IF ( z0h_water == 9999999.9_wp )  THEN
-             z0h_water = water_pars(2,water_type)    
-          ENDIF  
-
-          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
-             albedo_type = INT(water_pars(3,water_type))        
-          ENDIF
-   
-          IF ( emissivity == 9999999.9_wp )  THEN
-             emissivity = water_pars(4,water_type)        
-          ENDIF
-
-       ENDIF 
-       
-       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(0,pavement_type) 
-          ENDIF
-
-          IF ( z0h_pavement == 9999999.9_wp )  THEN
-             z0h_pavement = pavement_pars(1,pavement_type)
-          ENDIF
-
-          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
-             pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
-          ENDIF
-
-          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
-             pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
-          ENDIF   
-    
-          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
-             albedo_type = INT(pavement_pars(2,pavement_type))        
-          ENDIF
-   
-          IF ( emissivity == 9999999.9_wp )  THEN
-             emissivity = pavement_pars(3,pavement_type)        
-          ENDIF
-
-!
-!--       If the depth level of the pavement is not set, determine it from
-!--       lookup table.
-          IF ( pavement_depth_level == 0 )  THEN
-             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)  THEN
-                   nzt_pavement = k-1
-                   EXIT
-                ENDIF
-             ENDDO
-          ELSE
-             nzt_pavement = pavement_depth_level
           ENDIF
 
@@ -2494,76 +2628,1348 @@
           surf_lsm_v(l)%r_soil_min    = min_soil_resistance
        ENDDO
-
+!
+!--    Level 2, initialization of soil parameters via soil_type read from file.
+!--    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 
+!--       vertical dimension is assumed.
+          IF ( soil_type_f%lod == 1 )  THEN
+!
+!--          Horizontal surfaces
+             DO  m = 1, surf_lsm_h%ns
+                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
+                   surf_lsm_h%alpha_vg(:,m)    = soil_pars(0,st)
+                   surf_lsm_h%l_vg(:,m)        = soil_pars(1,st)
+                   surf_lsm_h%n_vg(:,m)        = soil_pars(2,st)
+                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st)
+                   surf_lsm_h%m_sat(:,m)       = soil_pars(4,st)
+                   surf_lsm_h%m_fc(:,m)        = soil_pars(5,st)
+                   surf_lsm_h%m_wilt(:,m)      = soil_pars(6,st)
+                   surf_lsm_h%m_res(:,m)       = soil_pars(7,st)
+                ENDIF
+             ENDDO
+!
+!--          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) ) 
+                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
+                                                   surf_lsm_v(l)%building_covered(m) ) 
+
+                   st = soil_type_f%var_2d(j,i)
+                   IF ( st /= soil_type_f%fill )  THEN
+                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars(0,st)
+                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars(1,st)
+                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars(2,st)
+                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st)
+                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars(4,st)
+                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars(5,st)
+                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars(6,st)
+                      surf_lsm_v(l)%m_res(:,m)       = soil_pars(7,st)
+                   ENDIF
+                ENDDO
+             ENDDO
+!
+!--       Level of detail = 2, i.e. soil type and thus the soil parameters
+!--       can be heterogeneous along the vertical dimension. 
+          ELSE
+!
+!--          Horizontal surfaces
+             DO  m = 1, surf_lsm_h%ns
+                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)
+                   IF ( st /= soil_type_f%fill )  THEN
+                      surf_lsm_h%alpha_vg(k,m)    = soil_pars(0,st)
+                      surf_lsm_h%l_vg(k,m)        = soil_pars(1,st)
+                      surf_lsm_h%n_vg(k,m)        = soil_pars(2,st)
+                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st)
+                      surf_lsm_h%m_sat(k,m)       = soil_pars(4,st)
+                      surf_lsm_h%m_fc(k,m)        = soil_pars(5,st)
+                      surf_lsm_h%m_wilt(k,m)      = soil_pars(6,st)
+                      surf_lsm_h%m_res(k,m)       = soil_pars(7,st)
+                   ENDIF
+                ENDDO
+             ENDDO
+!
+!--          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) ) 
+                   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
+                      st = soil_type_f%var_3d(k,j,i)
+                      IF ( st /= soil_type_f%fill )  THEN
+                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars(0,st)
+                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars(1,st)
+                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars(2,st)
+                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st)
+                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars(4,st)
+                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars(5,st)
+                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars(6,st)
+                         surf_lsm_v(l)%m_res(k,m)       = soil_pars(7,st)
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+       ENDIF
+!
+!--    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 
+!--       parameters is assumed.
+!--       Horizontal surfaces
+          IF ( soil_pars_f%lod == 1 )  THEN
+!
+!--          Horizontal surfaces 
+             DO  m = 1, surf_lsm_h%ns
+                i = surf_lsm_h%i(m)
+                j = surf_lsm_h%j(m)
+
+                IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
+                IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
+                IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
+                IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
+                IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
+                IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
+                IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
+                IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )              &
+                   surf_lsm_h%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
+
+             ENDDO
+!
+!--          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) ) 
+                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
+                                                   surf_lsm_v(l)%building_covered(m) ) 
+
+                   IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
+                   IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
+                   IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
+                   IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
+                   IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
+                   IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
+                   IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
+                   IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )           &
+                      surf_lsm_v(l)%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
+
+                ENDDO
+             ENDDO
+!
+!--       Level of detail = 2, i.e. soil parameters can be set at each soil 
+!--       layer individually. 
+          ELSE
+!
+!--          Horizontal surfaces
+             DO  m = 1, surf_lsm_h%ns
+                i = surf_lsm_h%i(m)
+                j = surf_lsm_h%j(m)
+
+                DO  k = nzb_soil, nzt_soil
+                   IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
+                   IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
+                   IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
+                   IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
+                   IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
+                   IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
+                   IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
+                   IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
+                      surf_lsm_h%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
+                ENDDO
+
+             ENDDO
+!
+!--          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) ) 
+                   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
+                      IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
+                      IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
+                      IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
+                      IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
+                      IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
+                      IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
+                      IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
+                      IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
+                         surf_lsm_v(l)%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
+                   ENDDO
+
+                ENDDO
+             ENDDO
+
+          ENDIF
+       ENDIF
+
+!
+!--    Level 1, initialization of vegetation parameters. A horizontally 
+!--    homogeneous distribution is assumed here. 
+       IF ( vegetation_type /= 0 )  THEN
+
+          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
+             min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type)
+          ENDIF
+
+          IF ( leaf_area_index == 9999999.9_wp )  THEN
+             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)     
+          ENDIF
+
+          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
+              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) 
+          ENDIF
+
+          IF ( z0h_vegetation == 9999999.9_wp )  THEN
+             z0h_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) 
+          ENDIF
+
+          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
+             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)        
+          ENDIF
+
+          IF ( c_surface == 9999999.9_wp )  THEN
+             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
+   
+          IF ( emissivity == 9999999.9_wp )  THEN
+             emissivity = vegetation_pars(ind_v_emis,vegetation_type)     
+          ENDIF
+
+       ENDIF
+!
+!--    Map values onto horizontal elemements
+       DO  m = 1, surf_lsm_h%ns
+          IF ( surf_lsm_h%vegetation_surface(m) )  THEN
+             surf_lsm_h%r_canopy_min(m)     = min_canopy_resistance
+             surf_lsm_h%lai(m)              = leaf_area_index
+             surf_lsm_h%c_veg(m)            = vegetation_coverage
+             surf_lsm_h%g_d(m)              = canopy_resistance_coefficient
+             surf_lsm_h%z0(m)               = z0_vegetation
+             surf_lsm_h%z0h(m)              = z0h_vegetation
+             surf_lsm_h%z0q(m)              = z0h_vegetation
+             surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
+             surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
+             surf_lsm_h%f_sw_in(m)          = f_shortwave_incoming
+             surf_lsm_h%c_surface(m)        = c_surface
+             surf_lsm_h%albedo_type(ind_veg,m) = albedo_type
+             surf_lsm_h%emissivity(ind_veg,m)  = emissivity
+          ELSE
+             surf_lsm_h%lai(m)   = 0.0_wp
+             surf_lsm_h%c_veg(m) = 0.0_wp
+             surf_lsm_h%g_d(m)   = 0.0_wp
+          ENDIF
+  
+       ENDDO
+!
+!--    Map values onto vertical elements, even though this does not make
+!--    much sense.
+       DO  l = 0, 3
+          DO  m = 1, surf_lsm_v(l)%ns
+             IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
+                surf_lsm_v(l)%r_canopy_min(m)     = min_canopy_resistance
+                surf_lsm_v(l)%lai(m)              = leaf_area_index
+                surf_lsm_v(l)%c_veg(m)            = vegetation_coverage
+                surf_lsm_v(l)%g_d(m)              = canopy_resistance_coefficient
+                surf_lsm_v(l)%z0(m)               = z0_vegetation
+                surf_lsm_v(l)%z0h(m)              = z0h_vegetation
+                surf_lsm_v(l)%z0q(m)              = z0h_vegetation
+                surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
+                surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
+                surf_lsm_v(l)%f_sw_in(m)          = f_shortwave_incoming
+                surf_lsm_v(l)%c_surface(m)        = c_surface
+                surf_lsm_v(l)%albedo_type(ind_veg,m) = albedo_type
+                surf_lsm_v(l)%emissivity(ind_veg,m)  = emissivity
+             ELSE
+                surf_lsm_v(l)%lai(m)   = 0.0_wp
+                surf_lsm_v(l)%c_veg(m) = 0.0_wp
+                surf_lsm_v(l)%g_d(m)   = 0.0_wp
+             ENDIF
+          ENDDO
+       ENDDO
+
+!
+!--    Level 2, initialization of vegation parameters via vegetation_type read
+!--    from file. Vegetation parameters are initialized for each (y,x)-grid point 
+!--    individually using default paramter settings according to the given
+!--    vegetation type.
+       IF ( vegetation_type_f%from_file )  THEN
+!
+!--       Horizontal surfaces
+          DO  m = 1, surf_lsm_h%ns
+             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
+                surf_lsm_h%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
+                surf_lsm_h%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
+                surf_lsm_h%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
+                surf_lsm_h%g_d(m)              = vegetation_pars(ind_v_gd,st)
+                surf_lsm_h%z0(m)               = vegetation_pars(ind_v_z0,st)
+                surf_lsm_h%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
+                surf_lsm_h%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
+                surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
+                surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
+                surf_lsm_h%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
+                surf_lsm_h%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
+                surf_lsm_h%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) )
+                surf_lsm_h%emissivity(ind_veg,m)  = vegetation_pars(ind_v_emis,st)
+             ENDIF
+          ENDDO
+!
+!--       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) ) 
+                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
+                                                   surf_lsm_v(l)%building_covered(m) ) 
+             
+                st = vegetation_type_f%var(j,i)
+                IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
+                   surf_lsm_v(l)%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
+                   surf_lsm_v(l)%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
+                   surf_lsm_v(l)%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
+                   surf_lsm_v(l)%g_d(m)              = vegetation_pars(ind_v_gd,st)
+                   surf_lsm_v(l)%z0(m)               = vegetation_pars(ind_v_z0,st)
+                   surf_lsm_v(l)%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
+                   surf_lsm_v(l)%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
+                   surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
+                   surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
+                   surf_lsm_v(l)%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
+                   surf_lsm_v(l)%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
+                   surf_lsm_v(l)%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) )
+                   surf_lsm_v(l)%emissivity(ind_veg,m)  = vegetation_pars(ind_v_emis,st)
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    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 
+          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 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.   
+             IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.                &
+                   ANY( vegetation_pars_f%pars_xy(:,j,i) /=                    &
+                   vegetation_pars_f%fill ) )  THEN
+                WRITE( message_string, * )                                     &
+                                 'surface element at grid point (j,i) = (',    &
+                                 j, i, ') is not a vegation surface, ',        &
+                                 'so that information given in ',              &
+                                 'vegetation_pars at this point is neglected.' 
+                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
+             ELSE
+
+                IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=            &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%r_canopy_min(m)  =                               &
+                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=            &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%lai(m)           =                               &
+                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=             &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%c_veg(m)         =                               &
+                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=                &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%g_d(m)           =                               &
+                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=                &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%z0(m)            =                               &
+                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=              &
+                     vegetation_pars_f%fill )  THEN
+                   surf_lsm_h%z0h(m)           =                               &
+                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
+                   surf_lsm_h%z0q(m)           =                               &
+                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
+                ENDIF
+                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=          &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%lambda_surface_s(m) =                            &
+                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=          &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%lambda_surface_u(m) =                            &
+                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=           &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%f_sw_in(m)          =                            &
+                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=            &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%c_surface(m)        =                            &
+                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
+                IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=                &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%albedo_type(ind_veg,m) =                         &
+                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
+                IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=              &
+                     vegetation_pars_f%fill )                                  &
+                   surf_lsm_h%emissivity(ind_veg,m)  =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
+             ENDIF
+          ENDDO
+!
+!--       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) ) 
+                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 
+!--             vegetation_pars is given, neglect this information and give an
+!--             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) /=                 &
+                      vegetation_pars_f%fill ) )  THEN
+                   WRITE( message_string, * )                                  &
+                                 'surface element at grid point (j,i) = (',    &
+                                 j, i, ') is not a vegation surface, ',        &
+                                 'so that information given in ',              &
+                                 'vegetation_pars at this point is neglected.' 
+                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
+                ELSE
+
+                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=         &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%r_canopy_min(m)  =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=         &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%lai(m)           =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=          &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%c_veg(m)         =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=             &
+                        vegetation_pars_f%fill )                               &
+                     surf_lsm_v(l)%g_d(m)            =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=             &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%z0(m)            =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=           &
+                        vegetation_pars_f%fill )  THEN
+                      surf_lsm_v(l)%z0h(m)           =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
+                      surf_lsm_v(l)%z0q(m)           =                         &
+                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
+                   ENDIF
+                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=       &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%lambda_surface_s(m)  =                     &
+                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=       &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%lambda_surface_u(m)  =                     &
+                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=        &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%f_sw_in(m)           =                     &
+                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=         &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%c_surface(m)         =                     &
+                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
+                   IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=             &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%albedo_type(ind_veg,m) =                   &
+                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
+                   IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=           &
+                        vegetation_pars_f%fill )                               &
+                      surf_lsm_v(l)%emissivity(ind_veg,m)  =                   &
+                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
+                ENDIF
+
+             ENDDO
+          ENDDO
+       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)       
+          ENDIF
+
+          IF ( z0_water == 9999999.9_wp )  THEN
+             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  
+
+          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
+             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 
+!
+!--    Map values onto horizontal elemements
+       DO  m = 1, surf_lsm_h%ns
+          IF ( surf_lsm_h%water_surface(m) )  THEN
+             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )        &
+                t_soil_h%var_2d(:,m)           = water_temperature
+             surf_lsm_h%z0(m)               = z0_water
+             surf_lsm_h%z0h(m)              = z0h_water
+             surf_lsm_h%z0q(m)              = z0h_water
+             surf_lsm_h%lambda_surface_s(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(ind_wat,m) = albedo_type
+             surf_lsm_h%emissivity(ind_wat,m)  = emissivity
+          ENDIF
+       ENDDO
+!
+!--    Map values onto vertical elements, even though this does not make
+!--    much sense.
+       DO  l = 0, 3
+          DO  m = 1, surf_lsm_v(l)%ns
+             IF ( surf_lsm_v(l)%water_surface(m) )  THEN
+                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
+                   t_soil_v(l)%var_2d(:,m)           = water_temperature
+                surf_lsm_v(l)%z0(m)               = z0_water
+                surf_lsm_v(l)%z0h(m)              = z0h_water
+                surf_lsm_v(l)%z0q(m)              = z0h_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)%c_surface(m)        = 0.0_wp
+                surf_lsm_v(l)%albedo_type(ind_wat,m) = albedo_type
+                surf_lsm_v(l)%emissivity(ind_wat,m)  = emissivity
+             ENDIF 
+          ENDDO
+       ENDDO
+!
+!
+!--    Level 2, initialization of water parameters via water_type read
+!--    from file. Water surfaces are initialized for each (y,x)-grid point 
+!--    individually using default paramter settings according to the given
+!--    water type.
+!--    Note, parameter 3/4 of water_pars are albedo and emissivity,
+!--    whereas paramter 3/4 of water_pars_f are heat conductivities!
+       IF ( water_type_f%from_file )  THEN
+!
+!--       Horizontal surfaces
+          DO  m = 1, surf_lsm_h%ns
+             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
+                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
+                   t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st)
+                surf_lsm_h%z0(m)     = water_pars(ind_w_z0,st)
+                surf_lsm_h%z0h(m)    = water_pars(ind_w_z0h,st)
+                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%c_surface(m)        = 0.0_wp
+                surf_lsm_h%albedo_type(ind_wat,m) = INT( water_pars(ind_w_at,st) )
+                surf_lsm_h%emissivity(ind_wat,m)  = water_pars(ind_w_emis,st)
+             ENDIF
+          ENDDO
+!
+!--       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) ) 
+                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
+                                                surf_lsm_v(l)%building_covered(m) ) 
+             
+                st = water_type_f%var(j,i)
+                IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
+                   IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  &
+                      t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st)
+                   surf_lsm_v(l)%z0(m)     = water_pars(ind_w_z0,st)
+                   surf_lsm_v(l)%z0h(m)    = water_pars(ind_w_z0h,st)
+                   surf_lsm_v(l)%z0q(m)    = water_pars(ind_w_z0h,st)
+                   surf_lsm_v(l)%lambda_surface_s(m) =                         &
+                                                   water_pars(ind_w_lambda_s,st)
+                   surf_lsm_v(l)%lambda_surface_u(m) =                         &
+                                                   water_pars(ind_w_lambda_u,st)           
+                   surf_lsm_v(l)%c_surface(m)     = 0.0_wp
+                   surf_lsm_v(l)%albedo_type(ind_wat,m) =                      &
+                                                  INT( water_pars(ind_w_at,st) )
+                   surf_lsm_v(l)%emissivity(ind_wat,m)  =                      &
+                                                  water_pars(ind_w_emis,st)
+                ENDIF
+             ENDDO
+          ENDDO
+       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 
+          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 
+!--          water_pars is given, neglect this information and give an
+!--          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
+                WRITE( message_string, * )                                     &
+                              'surface element at grid point (j,i) = (',       &
+                              j, i, ') is not a water surface, ',              &
+                              'so that information given in ',                 &
+                              'water_pars at this point is neglected.' 
+                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
+             ELSE
+                IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                   &
+                     water_pars_f%fill  .AND.                                  &
+                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
+                      t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
+
+                IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) &
+                   surf_lsm_h%z0(m)     = water_pars_f%pars_xy(ind_w_z0,j,i)
+
+                IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )&
+                THEN
+                   surf_lsm_h%z0h(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
+                   surf_lsm_h%z0q(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
+                ENDIF
+
+                IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=               &
+                     water_pars_f%fill )                                       &
+                   surf_lsm_h%lambda_surface_s(m) =                            &
+                                        water_pars_f%pars_xy(ind_w_lambda_s,j,i)
+
+                IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=               &
+                      water_pars_f%fill )                                      &
+                   surf_lsm_h%lambda_surface_u(m) =                            &
+                                        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 )                                       &
+                   surf_lsm_h%albedo_type(ind_wat,m) =                         &
+                                       INT( water_pars_f%pars_xy(ind_w_at,j,i) )
+
+                IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                   &
+                     water_pars_f%fill )                                       &
+                   surf_lsm_h%emissivity(ind_wat,m) =                          &
+                   water_pars_f%pars_xy(ind_w_emis,j,i)  
+             ENDIF
+          ENDDO
+!
+!--       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) ) 
+                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 
+!--             water_pars is given, neglect this information and give an
+!--             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) /=                      &
+                      water_pars_f%fill ) )  THEN
+                   WRITE( message_string, * )                                  &
+                              'surface element at grid point (j,i) = (',       &
+                              j, i, ') is not a water surface, ',              &
+                              'so that information given in ',                 &
+                              'water_pars at this point is neglected.' 
+                   CALL message( 'land_surface_model_mod', 'PA0999',           &
+                                  0, 0, 0, 6, 0 )
+                ELSE
+
+                   IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                &
+                     water_pars_f%fill  .AND.                                  &
+                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
+                      t_soil_v(l)%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
+
+                   IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /=                  &
+                        water_pars_f%fill )                                    &
+                      surf_lsm_v(l)%z0(m)   = water_pars_f%pars_xy(ind_w_z0,j,i)
+
+                   IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /=                 &
+                       water_pars_f%fill )  THEN
+                      surf_lsm_v(l)%z0h(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
+                      surf_lsm_v(l)%z0q(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
+                   ENDIF
+
+                   IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=            &
+                        water_pars_f%fill )                                    &
+                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
+                                      water_pars_f%pars_xy(ind_w_lambda_s,j,i)
+
+                   IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=            &
+                        water_pars_f%fill )                                    &
+                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
+                                      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 )                                    &
+                      surf_lsm_v(l)%albedo_type(ind_wat,m)    =                &
+                                      INT( water_pars_f%pars_xy(ind_w_at,j,i) )
+
+                   IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                &
+                        water_pars_f%fill )                                    &
+                      surf_lsm_v(l)%emissivity(ind_wat,m)     =                &
+                                      water_pars_f%pars_xy(ind_w_emis,j,i)  
+                ENDIF
+             ENDDO
+          ENDDO
+
+       ENDIF
+!
+!--    Initialize pavement-type surfaces, level 1
+       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) 
+          ENDIF
+
+          IF ( z0h_pavement == 9999999.9_wp )  THEN
+             z0h_pavement = pavement_pars(ind_p_z0h,pavement_type)
+          ENDIF
+
+          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
+             pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
+          ENDIF
+
+          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
+             pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
+          ENDIF   
+    
+          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
+             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)        
+          ENDIF
+
+!
+!--       If the depth level of the pavement is not set, determine it from
+!--       lookup table.
+          IF ( pavement_depth_level == 0 )  THEN
+             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)&
+                THEN
+                   nzt_pavement = k-1
+                   EXIT
+                ENDIF
+             ENDDO
+          ELSE
+             nzt_pavement = pavement_depth_level
+          ENDIF
+
+       ENDIF
+!
+!--    Level 1 initialization of pavement type surfaces. Horizontally 
+!--    homogeneous characteristics are assumed
+       surf_lsm_h%nzt_pavement = pavement_depth_level
+       DO  m = 1, surf_lsm_h%ns
+          IF ( surf_lsm_h%pavement_surface(m) )  THEN
+             surf_lsm_h%nzt_pavement(m)        = nzt_pavement
+             surf_lsm_h%z0(m)                  = z0_pavement
+             surf_lsm_h%z0h(m)                 = z0h_pavement
+             surf_lsm_h%z0q(m)                 = z0h_pavement
+             surf_lsm_h%lambda_surface_s(m)    = pavement_heat_conduct         &
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp    
+             surf_lsm_h%lambda_surface_u(m)    = pavement_heat_conduct         &
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp           
+             surf_lsm_h%c_surface(m)           = pavement_heat_capacity        &
+                                                        * dz_soil(nzb_soil)    &
+                                                        * 0.25_wp                                   
+
+             surf_lsm_h%albedo_type(ind_pav,m) = albedo_type
+             surf_lsm_h%emissivity(ind_pav,m)  = emissivity      
+      
+             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)                       
+                   surf_lsm_h%rho_c_total_def(k,m) =                           &
+                                     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
+       ENDDO                               
+
+       DO  l = 0, 3
+          surf_lsm_v(l)%nzt_pavement = pavement_depth_level
+          DO  m = 1, surf_lsm_v(l)%ns
+             IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
+                surf_lsm_v(l)%nzt_pavement(m)        = nzt_pavement
+                surf_lsm_v(l)%z0(m)                  = z0_pavement
+                surf_lsm_v(l)%z0h(m)                 = z0h_pavement
+                surf_lsm_v(l)%z0q(m)                 = z0h_pavement
+                surf_lsm_v(l)%lambda_surface_s(m)    = pavement_heat_conduct   &
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp    
+                surf_lsm_v(l)%lambda_surface_u(m)    = pavement_heat_conduct   &
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp           
+                surf_lsm_v(l)%c_surface(m)           = pavement_heat_capacity  &
+                                                        * dz_soil(nzb_soil)    &
+                                                        * 0.25_wp                                      
+
+                surf_lsm_v(l)%albedo_type(ind_pav,m) = albedo_type
+                surf_lsm_v(l)%emissivity(ind_pav,m)  = emissivity
+
+                IF ( pavement_type /= 0 )  THEN
+                   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)                       
+                      surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
+                                     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
+          ENDDO
+       ENDDO
+!
+!--    Level 2 initialization of pavement type surfaces via pavement_type read
+!--    from file. Pavement surfaces are initialized for each (y,x)-grid point 
+!--    individually.
+       IF ( pavement_type_f%from_file )  THEN
+!
+!--       Horizontal surfaces
+          DO  m = 1, surf_lsm_h%ns
+             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  
+                   IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp       &
+                   .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)      &
+                   THEN
+                      surf_lsm_h%nzt_pavement(m) = k-1
+                      EXIT
+                   ENDIF
+                ENDDO
+
+                surf_lsm_h%z0(m)                = pavement_pars(ind_p_z0,st)
+                surf_lsm_h%z0h(m)               = pavement_pars(ind_p_z0h,st)
+                surf_lsm_h%z0q(m)               = pavement_pars(ind_p_z0h,st)
+
+                surf_lsm_h%lambda_surface_s(m)  =                              &
+                                              pavement_subsurface_pars_1(0,st) &
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp    
+                surf_lsm_h%lambda_surface_u(m)  =                              &
+                                              pavement_subsurface_pars_1(0,st) &
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp        
+                surf_lsm_h%c_surface(m)         =                              &
+                                               pavement_subsurface_pars_2(0,st)&
+                                                        * dz_soil(nzb_soil)    &
+                                                        * 0.25_wp                               
+                surf_lsm_h%albedo_type(ind_pav,m) = INT( pavement_pars(ind_p_at,st) )
+                surf_lsm_h%emissivity(ind_pav,m)  = pavement_pars(ind_p_emis,st)  
+
+                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)                       
+                   surf_lsm_h%rho_c_total_def(k,m) =                           &
+                                     pavement_subsurface_pars_2(k,pavement_type) 
+                ENDDO   
+             ENDIF
+          ENDDO
+!
+!--       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) ) 
+                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
+                                                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  
+                      IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp    &
+                      .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)   &
+                      THEN
+                         surf_lsm_v(l)%nzt_pavement(m) = k-1
+                         EXIT
+                      ENDIF
+                   ENDDO
+
+                   surf_lsm_v(l)%z0(m)  = pavement_pars(ind_p_z0,st)
+                   surf_lsm_v(l)%z0h(m) = pavement_pars(ind_p_z0h,st)
+                   surf_lsm_v(l)%z0q(m) = pavement_pars(ind_p_z0h,st)
+
+                   surf_lsm_v(l)%lambda_surface_s(m)  =                        &
+                                              pavement_subsurface_pars_1(0,st) &
+                                                  * 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     
+
+                   surf_lsm_v(l)%c_surface(m)    =                             &
+                                           pavement_subsurface_pars_2(0,st)    &
+                                                        * dz_soil(nzb_soil)    &
+                                                        * 0.25_wp                                    
+                   surf_lsm_v(l)%albedo_type(ind_pav,m) =                      &
+                                              INT( pavement_pars(ind_p_at,st) )
+                   surf_lsm_v(l)%emissivity(ind_pav,m)  =                      &
+                                              pavement_pars(ind_p_emis,st)    
+
+                   DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
+                      surf_lsm_v(l)%lambda_h_def(k,m)    =                    &
+                                    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   
+                ENDIF
+             ENDDO
+          ENDDO
+       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 
+          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 
+!--          pavement_pars is given, neglect this information and give an
+!--          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) /=                      &
+                   pavement_pars_f%fill ) )  THEN
+                WRITE( message_string, * )                                     &
+                              'surface element at grid point (j,i) = (',       &
+                              j, i, ') is not a pavement surface, ',           &
+                              'so that information given in ',                 &
+                              'pavement_pars at this point is neglected.' 
+                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
+             ELSE
+                IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /=                  &
+                     pavement_pars_f%fill )                                    &
+                   surf_lsm_h%z0(m)  = pavement_pars_f%pars_xy(ind_p_z0,j,i)
+                IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /=                 &
+                     pavement_pars_f%fill )  THEN
+                   surf_lsm_h%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
+                   surf_lsm_h%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
+                ENDIF
+                IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i) &
+                     /= pavement_subsurface_pars_f%fill )  THEN
+                   surf_lsm_h%lambda_surface_s(m)  =                           &
+                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp
+                   surf_lsm_h%lambda_surface_u(m)  =                           &
+                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp    
+                ENDIF
+                IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) /= &
+                     pavement_subsurface_pars_f%fill )  THEN
+                   surf_lsm_h%c_surface(m)     =                               &
+                      pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)   &
+                                                  * dz_soil(nzb_soil)          &
+                                                  * 0.25_wp                                    
+                ENDIF
+                IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=                  &
+                     pavement_pars_f%fill )                                    &
+                   surf_lsm_h%albedo_type(ind_pav,m) =                         &
+                                              INT( pavement_pars(ind_p_at,st) )
+                IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /=                &
+                     pavement_pars_f%fill )                                    &
+                   surf_lsm_h%emissivity(ind_pav,m)  =                         &
+                                              pavement_pars(ind_p_emis,st) 
+             ENDIF 
+
+          ENDDO
+!
+!--       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) ) 
+                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 
+!--             pavement_pars is given, neglect this information and give an
+!--             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) /=                   &
+                      pavement_pars_f%fill ) )  THEN
+                   WRITE( message_string, * )                                  &
+                                 'surface element at grid point (j,i) = (',    &
+                                 j, i, ') is not a pavement surface, ',        &
+                                 'so that information given in ',              &
+                                 'pavement_pars at this point is neglected.' 
+                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
+                ELSE
+
+                   IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /=               &
+                        pavement_pars_f%fill )                                 &
+                      surf_lsm_v(l)%z0(m) = pavement_pars_f%pars_xy(ind_p_z0,j,i)
+                   IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /=              &
+                        pavement_pars_f%fill )  THEN
+                      surf_lsm_v(l)%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
+                      surf_lsm_v(l)%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
+                   ENDIF
+                   IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
+                        /= pavement_subsurface_pars_f%fill )  THEN
+                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
+                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp
+                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
+                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
+                                                  * ddz_soil(nzb_soil)         &
+                                                  * 2.0_wp    
+                   ENDIF
+                   IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) &
+                        /= pavement_subsurface_pars_f%fill )  THEN
+                      surf_lsm_v(l)%c_surface(m)    =                          &
+                         pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)&
+                                                  * dz_soil(nzb_soil)          &
+                                                  * 0.25_wp                                  
+                   ENDIF
+                   IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=               &
+                        pavement_pars_f%fill )                                 &
+                      surf_lsm_v(l)%albedo_type(ind_pav,m) =                   &
+                                            INT( pavement_pars(ind_p_at,st) )
+
+                   IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /=             &
+                        pavement_pars_f%fill )                                 &
+                      surf_lsm_v(l)%emissivity(ind_pav,m)  =                   &
+                                            pavement_pars(ind_p_emis,st)  
+                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. 
+       IF ( pavement_subsurface_pars_f%from_file )  THEN
+!
+!--       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
+
+                i = surf_lsm_h%i(m)
+                j = surf_lsm_h%j(m)
+
+                surf_lsm_h%nzt_pavement(m) = 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)
+                   surf_lsm_h%rho_c_total_def(k,m) =                           &
+                       pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
+                ENDDO
+
+             ENDIF
+          ENDDO
+          DO  l = 0, 3
+             DO  m = 1, surf_lsm_v(l)%ns
+                IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
+
+                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
+                                                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)%nzt_pavement(m) = 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)
+                      surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
+                       pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
+                   ENDDO
+
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDIF
 !
 !--    Initial run actions
        IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 !
-!--       First, for horizontal surfaces
-!
-!--       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
+!--       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
           DO  m = 1, surf_lsm_h%ns
-
-             t_soil_h%var_2d(:,m)    = 0.0_wp
-             m_soil_h%var_2d(:,m)    = 0.0_wp
-             m_liq_h%var_1d(m)       = 0.0_wp
-
-!--          Map user settings of T and q for each soil layer
-!--          (make sure that the soil moisture does not drop below the permanent
-!--          wilting point) -> problems with devision by zero)
-             IF ( surf_lsm_h%water_surface(m) )  THEN
-
-                surf_lsm_h%z0(m)          = z0_water
-                surf_lsm_h%z0h(m)         = z0h_water
-                surf_lsm_h%z0q(m)         = z0h_water
-                t_soil_h%var_2d(:,m)      = water_temperature
-                surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
-                surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
-                surf_lsm_h%c_surface(m)        = 0.0_wp
-                
-             ELSEIF ( surf_lsm_h%pavement_surface(m) )  THEN
-
-                surf_lsm_h%z0(m)          = z0_pavement
-                surf_lsm_h%z0h(m)         = z0h_pavement
-                surf_lsm_h%z0q(m)         = z0h_pavement   
+             IF ( surf_lsm_h%vegetation_surface(m)  .OR.                       &
+                  surf_lsm_h%pavement_surface(m) )  THEN
                 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)
+                   t_soil_h%var_2d(k,m) = soil_temperature(k)
+                   m_soil_h%var_2d(k,m) = soil_moisture(k)
                 ENDDO
                 t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
-                surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct         &
-                                                 * ddz_soil(nzb_soil)          &
-                                                 * 2.0_wp
-                surf_lsm_h%lambda_surface_u(m) = surf_lsm_h%lambda_surface_s(m)                                         
-                surf_lsm_h%c_surface(m)        = pavement_heat_capacity        &
-                                                 * dz_soil(nzb_soil)           &
-                                                 * 0.25_wp
-
-                IF ( pavement_type /= 0 )  THEN
+             ENDIF
+          ENDDO
+          DO  l = 0, 3
+             DO  m = 1, surf_lsm_v(l)%ns
+                IF ( surf_lsm_v(l)%vegetation_surface(m)  .OR.                 &
+                     surf_lsm_v(l)%pavement_surface(m) )  THEN
                    DO  k = nzb_soil, nzt_soil  
-                      surf_lsm_h%lambda_h_def(k,m)      =  pavement_subsurface_pars_1(k,pavement_type)                       
-                      surf_lsm_h%rho_c_total_def(k,m)   =  pavement_subsurface_pars_2(k,pavement_type) 
+                      t_soil_v(l)%var_2d(k,m) = soil_temperature(k)
+                      m_soil_v(l)%var_2d(k,m) = soil_moisture(k)
                    ENDDO
-                ELSE
-                    surf_lsm_h%lambda_h_def(:,m)     = pavement_heat_conduct
-                    surf_lsm_h%rho_c_total_def(:,m)  = pavement_heat_capacity
+                   t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
                 ENDIF
-
-             ELSEIF ( surf_lsm_h%vegetation_surface(m) )  THEN
-
-                surf_lsm_h%z0(m)          = z0_vegetation
-                surf_lsm_h%z0h(m)         = z0h_vegetation
-                surf_lsm_h%z0q(m)         = z0h_vegetation   
-                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)       
-                ENDDO 
-                t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
-                surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
-                surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
-                surf_lsm_h%c_surface(m)        = c_surface
-
+             ENDDO
+          ENDDO
+!
+!--       Level 2, if soil moisture and/or temperature  are 
+!--       provided from file, interpolate / extrapolate the provided data
+!--       onto the respective soil layers. Please note, both, zs as well as 
+!--       init_3d%z_soil indicate a depth with positive values, so that no 
+!--       distinction between atmosphere is required concerning interpolation.
+!--       Start with soil moisture
+          IF ( init_3d%from_file_msoil )  THEN
+
+             IF ( init_3d%lod_msoil == 1 )  THEN
+                DO  m = 1, surf_lsm_h%ns
+
+                   CALL netcdf_data_input_interpolate(                         &
+                                       m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
+                                       init_3d%msoil_init(:),                  &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                ENDDO
+                DO  l = 0, 3
+                   DO  m = 1, surf_lsm_v(l)%ns
+
+                      CALL netcdf_data_input_interpolate(                      &
+                                       m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
+                                       init_3d%msoil_init(:),                  &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                   ENDDO
+                ENDDO
+             ELSE
+
+                DO  m = 1, surf_lsm_h%ns
+                   i = surf_lsm_h%i(m)
+                   j = surf_lsm_h%j(m)
+
+                   IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )           &
+                      CALL netcdf_data_input_interpolate(                      &
+                                       m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
+                                       init_3d%msoil(:,j,i),                   &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                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) ) 
+                      j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,   &
+                                             surf_lsm_v(l)%building_covered(m) ) 
+
+                      IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )        &
+                         CALL netcdf_data_input_interpolate(                   &
+                                       m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
+                                       init_3d%msoil(:,j,i),                   &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                   ENDDO
+                ENDDO
              ENDIF
- 
+
+          ENDIF
+!
+!--       Soil temperature
+          IF ( init_3d%from_file_tsoil )  THEN
+
+             IF ( init_3d%lod_tsoil == 1 )  THEN ! change to 1 if provided correctly by INIFOR
+                DO  m = 1, surf_lsm_h%ns
+
+                   CALL netcdf_data_input_interpolate(                         &
+                                       t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
+                                       init_3d%tsoil_init(:),                  &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                   t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
+                ENDDO
+                DO  l = 0, 3
+                   DO  m = 1, surf_lsm_v(l)%ns
+
+                      CALL netcdf_data_input_interpolate(                      &
+                                       t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
+                                       init_3d%tsoil_init(:),                  &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                      t_soil_v(l)%var_2d(nzt_soil+1,m) =                       &
+                                                 t_soil_v(l)%var_2d(nzt_soil,m)
+                   ENDDO
+                ENDDO
+             ELSE
+
+                DO  m = 1, surf_lsm_h%ns
+                   i = surf_lsm_h%i(m)
+                   j = surf_lsm_h%j(m)
+
+                   IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )           &
+                      CALL netcdf_data_input_interpolate(                      &
+                                       t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
+                                       init_3d%tsoil(:,j,i),                   &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                   t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
+                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) ) 
+                      j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,   &
+                                                surf_lsm_v(l)%building_covered(m) ) 
+
+                      IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )        &
+                         CALL netcdf_data_input_interpolate(                   &
+                                       t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
+                                       init_3d%tsoil(:,j,i),                   &
+                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
+                                       nzb_soil, nzt_soil,                     &
+                                       nzb_soil, init_3d%nzs-1 )
+                      t_soil_v(l)%var_2d(nzt_soil+1,m) =                       &
+                                                 t_soil_v(l)%var_2d(nzt_soil,m)
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+!
+!--       Further initialization
+          DO  m = 1, surf_lsm_h%ns
+
              i   = surf_lsm_h%i(m)            
              j   = surf_lsm_h%j(m)
@@ -2575,92 +3981,35 @@
              IF ( surf_lsm_h%lambda_surface_s(m) == 0.0_wp )  THEN
                 t_surface_h%var_1d(:)    = t_soil_h%var_2d(nzb_soil,m)
-                surf_lsm_h%pt_surface(:) = t_soil_h%var_2d(nzb_soil,m) / exn
+                surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exn
              ELSE
-                t_surface_h%var_1d(:)    = pt_surface * exn
-                surf_lsm_h%pt_surface(:) = pt_surface  
+                t_surface_h%var_1d(m)    = pt(k-1,j,i) * exn
+                surf_lsm_h%pt_surface(m) = pt(k-1,j,i)  
              ENDIF
              
-             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
-                pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+             IF ( cloud_physics  .OR. cloud_droplets ) THEN
+                surf_lsm_h%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
              ELSE
-                pt1 = pt(k,j,i)
-             ENDIF
+                surf_lsm_h%pt1(m) = pt(k,j,i)
+             ENDIF 
+
 
 !
 !--          Assure that r_a cannot be zero at model start
-             IF ( pt1 == surf_lsm_h%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
+             IF ( surf_lsm_h%pt1(m) == surf_lsm_h%pt_surface(m) )              &
+                surf_lsm_h%pt1(m) = surf_lsm_h%pt1(m) + 1.0E-20_wp
 
              surf_lsm_h%us(m)   = 0.1_wp
-             surf_lsm_h%ts(m)   = ( pt1 - surf_lsm_h%pt_surface(m) )           &
+             surf_lsm_h%ts(m)   = ( surf_lsm_h%pt1(m) - surf_lsm_h%pt_surface(m) )&
                                   / surf_lsm_h%r_a(m)
              surf_lsm_h%shf(m)  = - surf_lsm_h%us(m) * surf_lsm_h%ts(m)        &
                                   * rho_surface
-
          ENDDO
-
 !
 !--      Vertical surfaces
          DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
-             
-                t_soil_v(l)%var_2d    = 0.0_wp
-                m_soil_v(l)%var_2d    = 0.0_wp
-                m_liq_v(l)%var_1d  = 0.0_wp
-
-             
-!--             Map user settings of T and q for each soil layer
-!--             (make sure that the soil moisture does not drop below the permanent
-!--             wilting point) -> problems with devision by zero)
-                IF ( surf_lsm_v(l)%water_surface(m) )  THEN
-                   surf_lsm_v(l)%z0(m)          = z0_water
-                   surf_lsm_v(l)%z0h(m)         = z0h_water 
-                   surf_lsm_v(l)%z0q(m)         = z0h_water                   
-                   t_soil_v(l)%var_2d(:,m)      = water_temperature
-                   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)%c_surface(m)        = 0.0_wp
-                   
-                ELSEIF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
-                   surf_lsm_v(l)%z0(m)          = z0_pavement
-                   surf_lsm_v(l)%z0h(m)         = z0h_pavement
-                   surf_lsm_v(l)%z0q(m)         = z0h_pavement  
-                   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) 
-                   ENDDO
-                   t_soil_v(l)%var_2d(nzt_soil+1,m)  = soil_temperature(nzt_soil+1)
-                   surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct   &
-                                                       * ddz_soil(nzb_soil)    &
-                                                       * 2.0_wp
-                   surf_lsm_v(l)%lambda_surface_u(m) = surf_lsm_v(l)%lambda_surface_s(m)
-                   surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity         &
-                                                * dz_soil(nzb_soil)            &
-                                                * 0.25_wp
-                IF ( pavement_type /= 0 )  THEN
-                   DO  k = nzb_soil, nzt_soil  
-                      surf_lsm_v(l)%lambda_h_def(k,m)    =  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
-                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
-                        
-                   
-                ELSEIF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
-
-                   surf_lsm_v(l)%z0(m)          = z0_vegetation
-                   surf_lsm_v(l)%z0h(m)         = z0h_vegetation
-                   surf_lsm_v(l)%z0q(m)         = z0h_vegetation   
-                   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)       
-                   ENDDO 
-                   t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
-                   surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
-                   surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
-                   surf_lsm_v(l)%c_surface(m)        = c_surface
-                ENDIF        
-         
+                i   = surf_lsm_v(l)%i(m)            
+                j   = surf_lsm_v(l)%j(m)
+                k   = surf_lsm_v(l)%k(m)         
 !
 !--             Calculate surface temperature. In case of bare soil, the surface 
@@ -2668,32 +4017,34 @@
 !--             layer
                 IF ( surf_lsm_v(l)%lambda_surface_s(m) == 0.0_wp )  THEN
-                   t_surface_v(l)%var_1d(:)      = t_soil_v(l)%var_2d(nzb_soil,m)
-                   surf_lsm_v(l)%pt_surface(:)   = t_soil_v(l)%var_2d(nzb_soil,m) / exn
+                   t_surface_v(l)%var_1d(m)      = t_soil_v(l)%var_2d(nzb_soil,m)
+                   surf_lsm_v(l)%pt_surface(m)   = t_soil_v(l)%var_2d(nzb_soil,m) / exn
                 ELSE
-                   t_surface_v(l)%var_1d(:)      = pt_surface * exn
-                   surf_lsm_v(l)%pt_surface(:)   = pt_surface               
+                   j_off = surf_lsm_v(l)%joff
+                   i_off = surf_lsm_v(l)%ioff
+
+                   t_surface_v(l)%var_1d(m)      = pt(k,j+j_off,i+i_off) * exn
+                   surf_lsm_v(l)%pt_surface(m)   = pt(k,j+j_off,i+i_off)           
                 ENDIF
 
+
+                IF ( cloud_physics  .OR. cloud_droplets ) THEN
+                   surf_lsm_v(l)%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+                ELSE
+                   surf_lsm_v(l)%pt1(m) = pt(k,j,i)
+                ENDIF 
+
+!
+!--             Assure that r_a cannot be zero at model start
+                IF ( surf_lsm_v(l)%pt1(m) == surf_lsm_v(l)%pt_surface(m) )     &
+                     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 
 !--             for the first time step. Only for interior core domain, not for ghost points
-
-                i   = surf_lsm_v(l)%i(m)            
-                j   = surf_lsm_v(l)%j(m)
-                k   = surf_lsm_v(l)%k(m)
-
-                IF ( cloud_physics  .OR.  cloud_droplets )  THEN
-                   pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
-                ELSE
-                   pt1 = pt(k,j,i)
-                ENDIF
-
-!
-!--             Assure that r_a cannot be zero at model start
-                IF ( pt1 == surf_lsm_v(l)%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
-
                 surf_lsm_v(l)%us(m)   = 0.1_wp
-                surf_lsm_v(l)%ts(m)   = ( pt1 - surf_lsm_v(l)%pt_surface(m) ) / surf_lsm_v(l)%r_a(m)
-                surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) * surf_lsm_v(l)%ts(m) * rho_surface
+                surf_lsm_v(l)%ts(m)   = ( surf_lsm_v(l)%pt1(m) - surf_lsm_v(l)%pt_surface(m) ) /&
+                                          surf_lsm_v(l)%r_a(m)
+                surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) *                &
+                                          surf_lsm_v(l)%ts(m) * rho_surface
+
              ENDDO
           ENDDO
@@ -2715,17 +4066,7 @@
 !
 !--          Set index offset of surface element, seen from atmospheric grid point 
-             IF ( l == 0 )  THEN 
-                j_off = -1 
-                i_off = 0
-             ELSEIF ( l == 1 )  THEN 
-                j_off = 1 
-                i_off = 0
-             ELSEIF ( l == 2 )  THEN 
-                j_off = 0 
-                i_off = -1
-             ELSEIF ( l == 3 )  THEN 
-                j_off = 0 
-                i_off = 1
-             ENDIF
+             j_off = surf_lsm_v(l)%joff
+             i_off = surf_lsm_v(l)%ioff
+
              DO  m = 1, surf_lsm_v(l)%ns
                 i   = surf_lsm_v(l)%i(m)            
@@ -2738,21 +4079,14 @@
        ENDIF
 !
-!--    Initialize root fraction
-!--    Horizontal surfaces
+!--    Level 1 initialization of root distribution - provided by the user via
+!--    via namelist.
        DO  m = 1, surf_lsm_h%ns
-          i   = surf_lsm_h%i(m)            
-          j   = surf_lsm_h%j(m)
-
           DO  k = nzb_soil, nzt_soil
              surf_lsm_h%root_fr(k,m) = root_fraction(k)
           ENDDO
        ENDDO
-!
-!--    Vertical surfaces
+
        DO  l = 0, 3
           DO  m = 1, surf_lsm_v(l)%ns
-             i   = surf_lsm_v(l)%i(m)            
-             j   = surf_lsm_v(l)%j(m)
-
              DO  k = nzb_soil, nzt_soil
                 surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
@@ -2762,8 +4096,8 @@
 
 !
+!--    Level 2 initialization of root distribution.
 !--    When no root distribution is given by the user, use look-up table to prescribe
-!--    the root fraction in the individual soil layers
+!--    the root fraction in the individual soil layers.
        IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
-
 !
 !--       First, calculate the index bounds for integration
@@ -2820,15 +4154,13 @@
 !--       Map calculated root fractions
           DO  m = 1, surf_lsm_h%ns
-             i   = surf_lsm_h%i(m)            
-             j   = surf_lsm_h%j(m)
              DO  k = nzb_soil, nzt_soil 
                 IF ( surf_lsm_h%pavement_surface(m)  .AND.                     &
-                     k <= nzt_pavement )  THEN
+                     k <= surf_lsm_h%nzt_pavement(m) )  THEN
                    surf_lsm_h%root_fr(k,m) = 0.0_wp
                 ELSE
                    surf_lsm_h%root_fr(k,m) = root_fraction(k)
                 ENDIF
-             ENDDO
-
+
+             ENDDO
 !
 !--          Normalize so that the sum = 1. Only relevant when the root         
@@ -2840,14 +4172,10 @@
                 ENDDO
              ENDIF
-
           ENDDO
           DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
-                i   = surf_lsm_v(l)%i(m)            
-                j   = surf_lsm_v(l)%j(m)
-
                 DO  k = nzb_soil, nzt_soil
                    IF ( surf_lsm_v(l)%pavement_surface(m)  .AND.               &
-                        k <= nzt_pavement )  THEN
+                        k <= surf_lsm_h%nzt_pavement(m) )  THEN
                       surf_lsm_v(l)%root_fr(k,m) = 0.0_wp
                    ELSE
@@ -2864,27 +4192,41 @@
                    ENDDO
                 ENDIF
-
              ENDDO
            ENDDO
        ENDIF
 !
-!--    Map vegetation parameters to the respective 2D arrays
-       surf_lsm_h%r_canopy_min      = min_canopy_resistance
-       surf_lsm_h%lai               = leaf_area_index
-       surf_lsm_h%c_veg             = vegetation_coverage
-       surf_lsm_h%g_d               = canopy_resistance_coefficient
-       surf_lsm_h%f_sw_in           = f_shortwave_incoming
-
-!--    Vertical surfaces
-       DO  l = 0, 3
-
-!
-!--       Map vegetation parameters to the respective 2D arrays
-          surf_lsm_v(l)%r_canopy_min      = min_canopy_resistance
-          surf_lsm_v(l)%lai               = leaf_area_index
-          surf_lsm_v(l)%c_veg             = vegetation_coverage
-          surf_lsm_v(l)%g_d               = canopy_resistance_coefficient
-          surf_lsm_v(l)%f_sw_in           = f_shortwave_incoming
-       ENDDO
+!--    Level 3 initialization of root distribution.
+!--    Take value from file
+       IF ( root_area_density_lsm_f%from_file )  THEN
+          DO  m = 1, surf_lsm_h%ns
+             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) ) 
+                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) 
+                ENDDO
+
+             ENDIF
+          ENDDO
+
+          DO  l = 0, 3
+             DO  m = 1, surf_lsm_v(l)%ns
+                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) ) 
+                   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) 
+                   ENDDO
+
+                ENDIF
+             ENDDO
+          ENDDO
+
+       ENDIF
  
 !
@@ -2912,8 +4254,32 @@
 !--    Store initial profiles of t_soil and m_soil (assuming they are 
 !--    horizontally homogeneous on this PE)
-       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),  & 
-                                                2, statistic_regions+1 )
+!--    DEACTIVATED FOR NOW - leads to error when number of locations with
+!--    soil model is zero on a PE.
+!        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),  & 
+!                                                 2, statistic_regions+1 )
+
+!
+!--    Finally, make some consistency checks. 
+!--    Ceck for eck for illegal combination of LAI and vegetation coverage.
+       IF ( ANY( .NOT. surf_lsm_h%pavement_surface  .AND.                      &
+                 surf_lsm_h%lai == 0.0_wp  .AND.  surf_lsm_h%c_veg == 1.0_wp ) &
+          )  THEN
+          message_string = 'For non-pavement surfaces the combination ' //     &
+                           ' lai = 0.0 and c_veg = 1.0 is not allowed.'
+          CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 )
+       ENDIF
+
+       DO  l = 0, 3
+          IF ( ANY( .NOT. surf_lsm_v(l)%pavement_surface  .AND.                &
+                    surf_lsm_v(l)%lai == 0.0_wp  .AND.                         &
+                    surf_lsm_v(l)%c_veg == 1.0_wp ) )  THEN
+             message_string = 'For non-pavement surfaces the combination ' //  &
+                              ' lai = 0.0 and c_veg = 1.0 is not allowed.'
+             CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 )
+          ENDIF
+       ENDDO
+
 
     END SUBROUTINE lsm_init
@@ -2943,11 +4309,8 @@
 !
 !--    Horizontal surfaces
-       ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns)                        )
        ALLOCATE ( m_liq_h_p%var_1d(1:surf_lsm_h%ns)                      )
        ALLOCATE ( t_surface_h%var_1d(1:surf_lsm_h%ns)                    )
        ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns)                  )
-       ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
        ALLOCATE ( m_soil_h_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
-       ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns)   )
        ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
 
@@ -2963,4 +4326,22 @@
           ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns)   )
           ALLOCATE ( t_soil_v_p(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
+       ENDDO
+!
+!--    Allocate soil temperature and moisture. As these variables might be 
+!--    already allocated in case of restarts, check this. 
+       IF ( .NOT. ALLOCATED( m_liq_h%var_1d ) )                                &
+          ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns) )
+       IF ( .NOT. ALLOCATED( m_soil_h%var_2d ) )                               &
+          ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
+       IF ( .NOT. ALLOCATED( t_soil_h%var_2d ) )                               &
+          ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
+
+       DO  l = 0, 3
+          IF ( .NOT. ALLOCATED( m_liq_v(l)%var_1d ) )                          &
+             ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns) )
+          IF ( .NOT. ALLOCATED( m_soil_v(l)%var_2d ) )                         &
+             ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
+          IF ( .NOT. ALLOCATED( t_soil_v(l)%var_2d ) )                         &
+             ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
        ENDDO
 #else
@@ -2988,4 +4369,12 @@
        ENDDO
 #endif
+!
+!--    Allocate array for heat flux in W/m2, required for radiation?
+!--    Consider to remove this array
+       ALLOCATE( surf_lsm_h%surfhf(1:surf_lsm_h%ns) )
+       DO  l = 0, 3
+          ALLOCATE( surf_lsm_v(l)%surfhf(1:surf_lsm_v(l)%ns) )
+       ENDDO
+
 
 !
@@ -3004,67 +4393,63 @@
           ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
        ENDDO 
-!
-!--    Allocate skin-surface temperature
-       ALLOCATE ( surf_lsm_h%pt_surface(1:surf_lsm_h%ns) )
-       DO  l = 0, 3
-          ALLOCATE ( surf_lsm_v(l)%pt_surface(1:surf_lsm_v(l)%ns) )
-       ENDDO 
 
 !
 !--    Allocate 2D vegetation model arrays
 !--    Horizontal surfaces
-       ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns) )
-       ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns)            )
-       ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns)        )
-       ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns)            )
-       ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns)          )
-       ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns)              )
-       ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns)              )
-       ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns)              )
-       ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns) )
-       ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns) )
+       ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns)    )
+       ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns)               )
+       ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns)           )
+       ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns)               )
+       ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns)             )
+       ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns)                 )
+       ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns)                 )
+       ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns)                 )
+       ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns)    )
+       ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns)    )
+       ALLOCATE ( surf_lsm_h%nzt_pavement(1:surf_lsm_h%ns)        )
+       ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns)    )
+       ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns)           )
+       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%r_a(1:surf_lsm_h%ns)                 )
+       ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)            )
+       ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns)              )
+       ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns)          )
+       ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)                 )
+       ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)        )
        ALLOCATE ( surf_lsm_h%vegetation_surface(1:surf_lsm_h%ns)  )
-       ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns) )
-       ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns)        )
-       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%r_a(1:surf_lsm_h%ns)              )
-       ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)         )
-       ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns)           )
-       ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns)       )
-       ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)              )
-       ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)     )
-       ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)    )
-
-       surf_lsm_h%water_surface     = .FALSE.
-       surf_lsm_h%pavement_surface  = .FALSE.
+       ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)       )
+
+       surf_lsm_h%water_surface        = .FALSE.
+       surf_lsm_h%pavement_surface     = .FALSE.
        surf_lsm_h%vegetation_surface   = .FALSE. 
 !
 !--    Vertical surfaces
        DO  l = 0, 3
-          ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns) )
-          ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns)            )
-          ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns)        )
-          ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns)            )
-          ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns)          )
-          ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns)              )
-          ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns)              )
-          ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns)              )
-          ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns) )
-          ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns) )
+          ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns)    )
+          ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns)               )
+          ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns)           )
+          ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns)               )
+          ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns)             )
+          ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns)                 )
+          ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns)                 )
+          ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns)                 )
+          ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns)    )
+          ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns)    )
+          ALLOCATE ( surf_lsm_v(l)%nzt_pavement(1:surf_lsm_v(l)%ns)        )
+          ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns)    )
+          ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns)           )
+          ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns)            )
+          ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns)            )
+          ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns)           )
+          ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)                 )
+          ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)            )
+          ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns)              )
+          ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns)          )
+          ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)                 )
+          ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)        )
           ALLOCATE ( surf_lsm_v(l)%vegetation_surface(1:surf_lsm_v(l)%ns)  )
-          ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns) )
-          ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns)        )
-          ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns)         )
-          ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns)         )
-          ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns)        )
-          ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)              )
-          ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)         )
-          ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns)           )
-          ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns)       )
-          ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)              )
-          ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)     )
-          ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)    )
+          ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)       )
 
           surf_lsm_v(l)%water_surface     = .FALSE.
@@ -3073,4 +4458,5 @@
           
        ENDDO
+
    
 #if ! defined( __nopointer )
@@ -3088,4 +4474,5 @@
        m_soil_v    => m_soil_v_1;    m_soil_v_p    => m_soil_v_2
        m_liq_v     => m_liq_v_1;     m_liq_v_p     => m_liq_v_2
+
 #endif
 
@@ -3122,5 +4509,6 @@
                                   residual_moisture, root_fraction,            &
                                   saturation_moisture, skip_time_do_lsm,       &
-                                  soil_moisture, soil_temperature, soil_type,  &
+                                  soil_moisture, soil_temperature,             &
+                                  soil_type,                                   &
                                   surface_type,                                &
                                   vegetation_coverage, vegetation_type,        &
@@ -3204,4 +4592,5 @@
        IF ( horizontal )  THEN
           surf           => surf_lsm_h
+
           surf_m_soil    => m_soil_h
           surf_m_soil_p  => m_soil_h_p
@@ -3212,4 +4601,5 @@
        ELSE
           surf           => surf_lsm_v(l)
+
           surf_m_soil    => m_soil_v(l)
           surf_m_soil_p  => m_soil_v_p(l)
@@ -3225,6 +4615,7 @@
              DO  k = nzb_soil, nzt_soil
 
-                IF ( surf%pavement_surface(m)  .AND.  k <= nzt_pavement )  THEN
-                    
+                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) 
@@ -3266,26 +4657,27 @@
              tend(:) = 0.0_wp
 
-             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *      &
-                            ( surf%lambda_h(nzb_soil,m)                        &
-                            * ( surf_t_soil%var_2d(nzb_soil+1,m)               &
-                            - surf_t_soil%var_2d(nzb_soil,m) )                 &
-                            * ddz_soil_center(nzb_soil) + surf%ghf(m) )        &
-                            * ddz_soil(nzb_soil)
+             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
+                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
+                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil) &
+                      + surf%ghf(m) ) * ddz_soil(nzb_soil)
 
              DO  k = nzb_soil+1, nzt_soil
-                tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )                      &
-                        * (   surf%lambda_h(k,m)                                  &
-                        * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
-                        * ddz_soil_center(k) - surf%lambda_h(k-1,m)               &
-                        * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
-                        * ddz_soil_center(k-1) ) * ddz_soil(k)
+                tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )                   &
+                          * (   surf%lambda_h(k,m)                             &
+                     * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
+                     * ddz_soil_center(k)                                      &
+                     - surf%lambda_h(k-1,m)                                    &
+                     * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
+                     * ddz_soil_center(k-1)                                    &
+                            ) * ddz_soil(k)
 
              ENDDO
 
              surf_t_soil_p%var_2d(nzb_soil:nzt_soil,m) =                       &
-                                  surf_t_soil%var_2d(nzb_soil:nzt_soil,m)      &
-                                  + dt_3d * ( tsc(2)                           &
-                                  * tend(nzb_soil:nzt_soil) + tsc(3)           &
-                                  * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
+                                       surf_t_soil%var_2d(nzb_soil:nzt_soil,m) &
+                                               + dt_3d * ( tsc(2)              &
+                                               * tend(nzb_soil:nzt_soil)       & 
+                                               + tsc(3)                        &
+                                               * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
 
 !
@@ -3299,7 +4691,7 @@
                          intermediate_timestep_count_max )  THEN
                    DO  k = nzb_soil, nzt_soil
-                      surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k)        &
-                                                   + 5.3125_wp                 &
-                                                   * surf_tt_soil_m%var_2d(k,m)
+                      surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) +      &
+                                                    5.3125_wp *                &
+                                                      surf_tt_soil_m%var_2d(k,m)
                    ENDDO
                 ENDIF
@@ -3312,11 +4704,11 @@
 !--             In order to prevent water tranport through paved surfaces, 
 !--             conductivity and diffusivity are set to zero
-                IF ( surf%pavement_surface(m)  .AND.  k <= nzt_pavement )  THEN
-
+                IF ( surf%pavement_surface(m)  .AND.                           &
+                     k <= surf%nzt_pavement(m) )  THEN
                    lambda_temp(k) = 0.0_wp
                    gamma_temp(k)  = 0.0_wp
-
+   
                 ELSE 
-
+    
 !
 !--                Calculate soil diffusivity at the center of the soil layers
@@ -3332,7 +4724,8 @@
                    h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) /          &
                               ( surf%m_res(k,m) -                              &
-                                MAX(surf_m_soil%var_2d(k,m), surf%m_wilt(k,m)) &
+                                MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )&
                               )                                                &
-                            )**( surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp )  &
+                            )**(                                               &
+                          surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp )         &
                                ) - 1.0_wp                                      &
                           )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m)
@@ -3340,9 +4733,13 @@
                    gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp +      &
                           ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m)        &
-                          )**( 1.0_wp - 1.0_wp / surf%n_vg(k,m) )              &
-                          - ( surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m)    &
-                          - 1.0_wp) )**2 ) / ( ( 1.0_wp + ( surf%alpha_vg(k,m) &
-                          * h_vg )**surf%n_vg(k,m) )**(( 1.0_wp  - 1.0_wp      &
-                              / surf%n_vg(k,m) ) * ( surf%l_vg(k,m) + 2.0_wp) ))
+                                                                  )**(         &
+                              1.0_wp - 1.0_wp / surf%n_vg(k,m)) - (            &
+                          surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m)        &
+                              - 1.0_wp) )**2 )                                 &
+                              / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg       &
+                              )**surf%n_vg(k,m) )**( ( 1.0_wp  - 1.0_wp        &
+                              / surf%n_vg(k,m) ) *                             &
+                              ( surf%l_vg(k,m) + 2.0_wp) ) )
+
                 ENDIF
 
@@ -3381,5 +4778,7 @@
 !--             would cause accumulation of (non-existing) water in the lowest
 !--             soil layer
-                IF ( conserve_water_content .AND. surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp )  THEN
+                IF ( conserve_water_content .AND.                              &
+                     surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp )  THEN
+
                    surf%gamma_w(nzt_soil,m) = 0.0_wp
                 ELSE
@@ -3404,5 +4803,5 @@
                     IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
                        m_total = m_total + surf%root_fr(k,m)                   &
-                                 * surf_m_soil%var_2d(k,m)
+                                         * surf_m_soil%var_2d(k,m)
                     ENDIF
                 ENDDO  
@@ -3411,5 +4810,5 @@
                       IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
                          root_extr(k) = surf%root_fr(k,m)                      &
-                                        * surf_m_soil%var_2d(k,m) / m_total
+                                      * surf_m_soil%var_2d(k,m) / m_total
                       ELSE
                          root_extr(k) = 0.0_wp
@@ -3428,4 +4827,5 @@
                             + surf%qsws_soil(m) ) * drho_l_lv )                &
                             * ddz_soil(nzb_soil)
+
 
                 DO  k = nzb_soil+1, nzt_soil-1
@@ -3451,6 +4851,7 @@
                 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) )   
+                                         + dt_3d * ( tsc(2) * tend(:)          &
+                                         + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
+   
 !
 !--             Account for dry soils (find a better solution here!)
@@ -3473,4 +4874,5 @@
                                                     * surf_tm_soil_m%var_2d(k,m)
                       ENDDO
+
                    ENDIF
                 ENDIF
@@ -3526,4 +4928,5 @@
                 m_liq_h   => m_liq_h_1;     m_liq_h_p      => m_liq_h_2
              ENDIF
+
 !
 !--          Vertical surfaces
@@ -3533,4 +4936,5 @@
                 m_soil_v  => m_soil_v_1;    m_soil_v_p     => m_soil_v_2
                 m_liq_v   => m_liq_v_1;     m_liq_v_p      => m_liq_v_2
+
              ENDIF
 
@@ -3545,4 +4949,5 @@
                 m_soil_h  => m_soil_h_2;    m_soil_h_p     => m_soil_h_1
                 m_liq_h   => m_liq_h_2;     m_liq_h_p      => m_liq_h_1
+
              ENDIF
 !
@@ -3891,4 +5296,5 @@
                 ENDDO
              ENDDO
+
 !
 !-- 
@@ -4327,5 +5733,5 @@
 !> Write restart data for land surface model
 !------------------------------------------------------------------------------!
- SUBROUTINE lsm_last_actions
+ SUBROUTINE lsm_write_restart_data
  
 
@@ -4335,4 +5741,7 @@
 
     IMPLICIT NONE
+
+    CHARACTER (LEN=1) ::  dum    !< dummy to create correct string for creating variable string
+    INTEGER(iwp)      ::  l      !< index variable for surface orientation
 
     IF ( write_binary )  THEN
@@ -4352,9 +5761,7 @@
           WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
        ENDIF
-       WRITE ( 14 )     'm_liq               ';  WRITE ( 14 )  m_liq_h%var_1d
        IF ( ALLOCATED( m_liq_av ) )  THEN
           WRITE ( 14 )  'm_liq_av            ';  WRITE ( 14 )  m_liq_av
        ENDIF
-       WRITE ( 14 )     'm_soil              ';  WRITE ( 14 )  m_soil_h%var_2d
        IF ( ALLOCATED( m_soil_av ) )  THEN
           WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
@@ -4369,14 +5776,71 @@
           WRITE ( 14 )  'qsws_veg_av         ';  WRITE ( 14 )  qsws_veg_av
        ENDIF 
-       WRITE ( 14 )     't_soil              ';  WRITE ( 14 )  t_soil_h%var_2d
        IF ( ALLOCATED( t_soil_av ) )  THEN
           WRITE ( 14 )  't_soil_av           ';  WRITE ( 14 )  t_soil_av
        ENDIF
 
+
+       WRITE ( 14 ) 'ns_h_on_file_lsm    '
+       WRITE ( 14 ) surf_lsm_h%ns
+       WRITE ( 14 ) 'ns_v_on_file_lsm    '
+       WRITE ( 14 ) surf_lsm_v(0:3)%ns
+       
+       WRITE ( 14 ) 'lsm_start_index_h   '
+       WRITE ( 14 ) surf_lsm_h%start_index
+       WRITE ( 14 ) 'lsm_end_index_h     '
+       WRITE ( 14 ) surf_lsm_h%end_index
+       WRITE ( 14 ) 't_soil_h            '
+       WRITE ( 14 ) t_soil_h%var_2d
+       
+       DO  l = 0, 3
+          WRITE ( 14 ) 'lsm_start_index_v   '
+          WRITE ( 14 ) surf_lsm_v(l)%start_index
+          WRITE ( 14 ) 'lsm_end_index_v     '
+          WRITE ( 14 ) surf_lsm_v(l)%end_index
+          WRITE( dum, '(I1)')  l          
+          WRITE ( 14 ) 't_soil_v(' // dum // ')         ' 
+          WRITE ( 14 ) t_soil_v(l)%var_2d          
+       ENDDO
+
+       WRITE ( 14 ) 'lsm_start_index_h   '
+       WRITE ( 14 ) surf_lsm_h%start_index
+       WRITE ( 14 ) 'lsm_end_index_h     '
+       WRITE ( 14 ) surf_lsm_h%end_index
+       WRITE ( 14 ) 'm_soil_h            '
+       WRITE ( 14 ) m_soil_h%var_2d
+       
+       DO  l = 0, 3
+          WRITE ( 14 ) 'lsm_start_index_v   '
+          WRITE ( 14 ) surf_lsm_v(l)%start_index
+          WRITE ( 14 ) 'lsm_end_index_v     '
+          WRITE ( 14 ) surf_lsm_v(l)%end_index
+          WRITE( dum, '(I1)')  l          
+          WRITE ( 14 ) 'm_soil_v(' // dum // ')         ' 
+          WRITE ( 14 ) m_soil_v(l)%var_2d          
+       ENDDO
+
+
+       WRITE ( 14 ) 'lsm_start_index_h   '
+       WRITE ( 14 ) surf_lsm_h%start_index
+       WRITE ( 14 ) 'lsm_end_index_h     '
+       WRITE ( 14 ) surf_lsm_h%end_index
+       WRITE ( 14 ) 'm_liq_h             '
+       WRITE ( 14 ) m_liq_h%var_1d
+       
+       DO  l = 0, 3
+          WRITE ( 14 ) 'lsm_start_index_v   '
+          WRITE ( 14 ) surf_lsm_v(l)%start_index
+          WRITE ( 14 ) 'lsm_end_index_v     '
+          WRITE ( 14 ) surf_lsm_v(l)%end_index
+          WRITE( dum, '(I1)')  l          
+          WRITE ( 14 ) 'm_liq_v(' // dum // ')          ' 
+          WRITE ( 14 ) m_liq_v(l)%var_1d          
+       ENDDO
+
        WRITE ( 14 )  '*** end lsm ***     '
 
     ENDIF
 
- END SUBROUTINE lsm_last_actions
+ END SUBROUTINE lsm_write_restart_data
 
 
@@ -4399,19 +5863,23 @@
     CHARACTER (LEN=20) :: field_char   !< 
 
-    INTEGER(iwp) ::  i               !< 
-    INTEGER(iwp) ::  k               !< 
-    INTEGER(iwp) ::  nxlc            !< 
-    INTEGER(iwp) ::  nxlf            !< 
-    INTEGER(iwp) ::  nxl_on_file     !< 
-    INTEGER(iwp) ::  nxrc            !< 
-    INTEGER(iwp) ::  nxrf            !< 
-    INTEGER(iwp) ::  nxr_on_file     !< 
-    INTEGER(iwp) ::  nync            !< 
-    INTEGER(iwp) ::  nynf            !< 
-    INTEGER(iwp) ::  nyn_on_file     !< 
-    INTEGER(iwp) ::  nysc            !< 
-    INTEGER(iwp) ::  nysf            !< 
-    INTEGER(iwp) ::  nys_on_file     !< 
-    INTEGER(iwp) ::  overlap_count   !< 
+    INTEGER(iwp) ::  i                 !< 
+    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) ::  nxl_on_file       !< 
+    INTEGER(iwp) ::  nxrc              !< 
+    INTEGER(iwp) ::  nxrf              !< 
+    INTEGER(iwp) ::  nxr_on_file       !< 
+    INTEGER(iwp) ::  nync              !< 
+    INTEGER(iwp) ::  nynf              !< 
+    INTEGER(iwp) ::  nyn_on_file       !< 
+    INTEGER(iwp) ::  nysc              !< 
+    INTEGER(iwp) ::  nysf              !< 
+    INTEGER(iwp) ::  nys_on_file       !< 
+    INTEGER(iwp) ::  overlap_count     !< 
+
+    INTEGER(iwp) ::  ns_v_on_file_lsm(0:3) !< number of vertical surface elements (natural type) on file
 
     INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !< 
@@ -4422,4 +5890,7 @@
     INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !< 
 
+    INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) ::  start_index_on_file
+    INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) ::  end_index_on_file
+
     REAL(wp),                                                                  &
        DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: &
@@ -4434,18 +5905,42 @@
           tmp_3d2   !< 
 
-    REAL(wp),                                                                  &
-       DIMENSION(1:surf_lsm_h%ns) ::                                           &
-          tmp_walltype_1d   !< 
-
-    REAL(wp),                                                                  &
-       DIMENSION(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ::                       &
-          tmp_walltype_2d   !< 
-
-    REAL(wp),                                                                  &
-       DIMENSION(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ::                         &
-          tmp_walltype_2d2  !< 
-
+    TYPE(surf_type_lsm) :: tmp_walltype_h_1d   !< temporary 1D array containing the respective surface variable stored on file, horizontal surfaces
+    TYPE(surf_type_lsm) :: tmp_walltype_h_2d   !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces
+    TYPE(surf_type_lsm) :: tmp_walltype_h_2d2  !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces 
+
+    TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_1d   !< temporary 1D array containing the respective surface variable stored on file, vertical surfaces 
+    TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_2d   !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces 
+    TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_2d2  !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces
 
    IF ( initializing_actions == 'read_restart_data' )  THEN
+      READ ( 13 )  field_char
+
+!
+!--   At first, determine the number of surface elements (with certain orientation) on file
+      IF ( TRIM( field_char ) /= 'ns_h_on_file_lsm' )  THEN
+!
+!--      Add a proper error message
+      ENDIF
+      READ ( 13 ) ns_h_on_file_lsm
+
+      READ ( 13 )  field_char
+      IF ( TRIM( field_char ) /= 'ns_v_on_file_lsm' )  THEN
+!
+!--      Add a proper error message
+      ENDIF
+      READ ( 13 ) ns_v_on_file_lsm
+
+!
+!--   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) )
+      ALLOCATE( tmp_walltype_h_2d2%var_2d(nzb_soil:nzt_soil,1:ns_h_on_file_lsm)  )
+
+      DO  l = 0, 3
+         ALLOCATE( tmp_walltype_v_1d(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)) )
+         ALLOCATE( tmp_walltype_v_2d2(l)%var_2d(nzb_soil:nzt_soil,1:ns_v_on_file_lsm(l))  )
+      ENDDO
+     
       READ ( 13 )  field_char
 
@@ -4463,4 +5958,5 @@
             nync = nynfa(i,k) + offset_ya(i,k)
 
+
             SELECT CASE ( TRIM( field_char ) )
 
@@ -4495,11 +5991,6 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   ghf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
+                   ghf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'm_liq' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_1d !tmp_2d
-                   m_liq_h%var_1d(1:surf_lsm_h%ns)  =                       &
-                                 tmp_walltype_1d(1:surf_lsm_h%ns)
 
                 CASE ( 'lai_av' )
@@ -4518,9 +6009,4 @@
                    m_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'm_soil' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d2(:,:)
-                   m_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
-                          tmp_walltype_2d2(:,1:surf_lsm_h%ns)   
 
                 CASE ( 'm_soil_av' )
@@ -4556,9 +6042,4 @@
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
-                CASE ( 't_soil' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d(:,:)
-                   t_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
-                         tmp_walltype_2d(:,1:surf_lsm_h%ns)      
-
                 CASE ( 't_soil_av' )
                    IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
@@ -4569,4 +6050,288 @@
                                     tmp_3d2(:,nysf-nbgp:nynf+nbgp,             &
                                     nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'lsm_start_index_h', 'lsm_start_index_v'  )   
+                   IF ( k == 1 )                                               &
+                      READ ( 13 )  start_index_on_file
+                      
+                CASE ( 'lsm_end_index_h', 'lsm_end_index_v' )   
+                   IF ( k == 1 )                                               &
+                      READ ( 13 )  end_index_on_file
+                
+                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
+                   CALL surface_restore_elements(                              &
+                                              t_soil_h%var_2d,                 &
+                                              tmp_walltype_h_2d%var_2d,        &
+                                              surf_lsm_h%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_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
+                   CALL surface_restore_elements(                              &
+                                           t_soil_v(0)%var_2d,                 &
+                                           tmp_walltype_v_2d(0)%var_2d,        &
+                                           surf_lsm_v(0)%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_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
+                   CALL surface_restore_elements(                              &
+                                           t_soil_v(1)%var_2d,                 &
+                                           tmp_walltype_v_2d(1)%var_2d,        &
+                                           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_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
+                   CALL surface_restore_elements(                              &
+                                           t_soil_v(2)%var_2d,                 &
+                                           tmp_walltype_v_2d(2)%var_2d,        &
+                                           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_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
+                   CALL surface_restore_elements(                              &
+                                           t_soil_v(3)%var_2d,                 &
+                                           tmp_walltype_v_2d(3)%var_2d,        &
+                                           surf_lsm_v(3)%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 ( '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
+                   CALL surface_restore_elements(                              &
+                                             m_soil_h%var_2d,                  &
+                                             tmp_walltype_h_2d2%var_2d,        &
+                                             surf_lsm_h%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 ( '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
+                   CALL surface_restore_elements(                              &
+                                          m_soil_v(0)%var_2d,                  & 
+                                          tmp_walltype_v_2d2(0)%var_2d,        &
+                                          surf_lsm_v(0)%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 ( '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
+                   CALL surface_restore_elements(                              &
+                                          m_soil_v(1)%var_2d,                  & 
+                                          tmp_walltype_v_2d2(1)%var_2d,        &
+                                          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 ( '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
+                   CALL surface_restore_elements(                              &
+                                          m_soil_v(2)%var_2d,                  & 
+                                          tmp_walltype_v_2d2(2)%var_2d,        &
+                                          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 ( '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
+                   CALL surface_restore_elements(                              &
+                                          m_soil_v(3)%var_2d,                  & 
+                                          tmp_walltype_v_2d2(3)%var_2d,        &
+                                          surf_lsm_v(3)%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 ( '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
+                   CALL surface_restore_elements(                              &
+                                              m_liq_h%var_1d,                  &
+                                              tmp_walltype_h_1d%var_1d,        &
+                                              surf_lsm_h%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 ( '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
+                   CALL surface_restore_elements(                              &
+                                           m_liq_v(0)%var_1d,                  &
+                                           tmp_walltype_v_1d(0)%var_1d,        &
+                                           surf_lsm_v(0)%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 ( '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
+                   CALL surface_restore_elements(                              &
+                                           m_liq_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 ( '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
+                   CALL surface_restore_elements(                              &
+                                           m_liq_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 ( '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
+                   CALL surface_restore_elements(                              &
+                                           m_liq_v(3)%var_1d,                  &
+                                           tmp_walltype_v_1d(3)%var_1d,        &
+                                           surf_lsm_v(3)%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 )
 
 
@@ -4645,5 +6410,5 @@
 !--          Charnock relation is used.
              surf_lsm_h%z0(m)  = ( 0.11_wp * molecular_viscosity /             &
-                                 surf_lsm_h%us(m) )  &
+                                 surf_lsm_h%us(m) )                            &
                                + ( alpha_ch * surf_lsm_h%us(m)**2 / g )
 
