Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2474)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2475)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: setting of vegetation_pars for bare soil corrected.
+! 
+! 2354 2017-08-17 10:49:36Z schwenkel
 ! minor bugfixes
 ! 
@@ -278,5 +281,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  albedo, albedo_type, emissivity, force_radiation_call, rad_net,         &
+        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
@@ -399,9 +402,9 @@
 
                                             
-    REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp,               & !< (NAMELIST) distribution of root surface area to the individual soil layers
-                                   soil_moisture = 0.0_wp,                     & !< NAMELIST soil moisture content (m3/m3)
-                                   soil_temperature = 300.0_wp,                & !< NAMELIST soil temperature (K) +1
-                                   dz_soil  = 9999999.9_wp,                    & !< (NAMELIST) soil layer depths (spacing)
-                                   zs_layer = 9999999.9_wp                       !< soil layer depths (edge)
+    REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp, & !< (NAMELIST) distribution of root surface area to the individual soil layers
+                                   soil_moisture = 0.0_wp,       & !< NAMELIST soil moisture content (m3/m3)
+                                   soil_temperature = 300.0_wp,  & !< NAMELIST soil temperature (K) +1
+                                   dz_soil  = 9999999.9_wp,      & !< (NAMELIST) soil layer depths (spacing)
+                                   zs_layer = 9999999.9_wp         !< soil layer depths (edge)
                                  
 #if defined( __nopointer )
@@ -412,8 +415,8 @@
 
     TYPE(surf_type_lsm), DIMENSION(0:3), TARGET  ::  &
-                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
-                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
-                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
-                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements
+                                     t_soil_v,       & !< Soil temperature (K), vertical surface elements
+                                     t_soil_v_p,     & !< Prog. soil temperature (K), vertical surface elements
+                                     m_soil_v,       & !< Soil moisture (m3/m3), vertical surface elements
+                                     m_soil_v_p        !< Prog. soil moisture (m3/m3), vertical surface elements
 #else
     TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
@@ -449,28 +452,28 @@
                                       t_surface_v,    & !< surface temperature (K), vertical surface elements 
                                       t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements 
-                                      m_liq_v,     & !< liquid water reservoir (m), vertical surface elements 
-                                      m_liq_v_p      !< progn. liquid water reservoir (m), vertical surface elements 
+                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements 
+                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements 
 #else
     TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements 
                                       t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements 
-                                      m_liq_h,     & !< liquid water reservoir (m), horizontal surface elements 
-                                      m_liq_h_p      !< progn. liquid water reservoir (m), horizontal surface elements 
+                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements 
+                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements 
 
     TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
                                       t_surface_h_2,  & !<
-                                      m_liq_h_1,   & !<
-                                      m_liq_h_2      !<
+                                      m_liq_h_1,      & !<
+                                      m_liq_h_2         !<
 
     TYPE(surf_type_lsm), DIMENSION(:), POINTER  ::    &
                                       t_surface_v,    & !< surface temperature (K), vertical surface elements 
                                       t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements 
-                                      m_liq_v,     & !< liquid water reservoir (m), vertical surface elements 
-                                      m_liq_v_p      !< progn. liquid water reservoir (m), vertical surface elements 
+                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements 
+                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements 
 
     TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
                                       t_surface_v_1,  & !<
                                       t_surface_v_2,  & !<
-                                      m_liq_v_1,   & !<
-                                      m_liq_v_2      !<
+                                      m_liq_v_1,      & !<
+                                      m_liq_v_2         !<
 #endif
 
@@ -489,10 +492,10 @@
 #endif
 
-    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m   !< liquid water reservoir tendency (m), horizontal surface elements 
+    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m      !< liquid water reservoir tendency (m), horizontal surface elements 
     TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements 
     TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements 
     TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements 
 
-    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m   !< liquid water reservoir tendency (m), vertical surface elements 
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m      !< liquid water reservoir tendency (m), vertical surface elements 
     TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements 
     TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements 
@@ -502,14 +505,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
                    
 
@@ -517,23 +520,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
                                                                  /)
 
@@ -541,11 +544,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
                                                                  /)
 
@@ -553,12 +556,12 @@
 !-- Pavement classes
     CHARACTER(20), DIMENSION(0:7), PARAMETER :: pavement_type_name = (/ &
-                                   'user defined        ',          & ! 0 
-                                   'asphalt             ',          & ! 1
-                                   'concrete            ',          & ! 2
-                                   'asphalt/concrete mix',          & ! 3
-                                   'brick pavers        ',          & ! 4
-                                   'cobblestone pavers  ',          & ! 5
-                                   'sett pavers         ',          & ! 6
-                                   'gravel pavers       '           & ! 7
+                                   'user defined        ', & ! 0 
+                                   'asphalt             ', & ! 1
+                                   'concrete            ', & ! 2
+                                   'asphalt/concrete mix', & ! 3
+                                   'brick pavers        ', & ! 4
+                                   'cobblestone pavers  ', & ! 5
+                                   'sett pavers         ', & ! 6
+                                   'gravel pavers       '  & ! 7
                                                                  /)                                                                 
                                                                  
@@ -566,10 +569,10 @@
 !-- 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
                                                                   /)                                                                  
                                                                  
@@ -581,5 +584,5 @@
 !-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface, albedo_type, emissivity
     REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
-          0.0_wp, 0.00_wp, 1.00_wp, 0.00_wp,  0.005_wp,   0.5E-4_wp,     0.0_wp,    0.0_wp, 0.00_wp, 0.00_wp,  0.0_wp, 0.94_wp, & !  1
+          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  0.005_wp,   0.5E-4_wp,     0.0_wp,    0.0_wp, 0.00_wp, 0.00_wp, 18.0_wp, 0.94_wp, & !  1
         180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp,   0.10_wp,    0.001_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.95_wp, & !  2
         110.0_wp, 2.00_wp, 1.00_wp, 0.00_wp,   0.03_wp,   0.3E-4_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.95_wp, & !  3
@@ -606,22 +609,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 /) )
 
@@ -630,13 +633,13 @@
 
 !
-!-- Soil parameters             alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res 
+!-- Soil parameters  alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res 
     REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/     &
-                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp,& ! 1
-                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp,& ! 2
-                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp,& ! 3
-                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp,& ! 4
-                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp,& ! 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 /) )
+                      3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp,& ! 1
+                      3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp,& ! 2
+                      0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp,& ! 3
+                      3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp,& ! 4
+                      2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp,& ! 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 /) )
 
 !
@@ -644,5 +647,5 @@
 !-- Pavement parameters  depth,        z0,       z0h, lambda_h,      rho_c, albedo_type, emissivity  
     REAL(wp), DIMENSION(0:6,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ &
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 17.0_wp, 0.97_wp, & ! 1
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 17.0_wp, 0.97_wp,  & ! 1
                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 18.0_wp, 0.94_wp,  & ! 2
                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 19.0_wp, 0.98_wp,  & ! 3                                  
@@ -650,16 +653,16 @@
                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 21.0_wp, 0.97_wp,  & ! 5
                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 22.0_wp, 0.97_wp,  & ! 6
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 23.0_wp, 0.97_wp  & ! 7
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 23.0_wp, 0.97_wp   & ! 7
                                                                  /), (/ 7, 7 /) )                             
 
 !
 !-- TO BE FILLED
-!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
+!-- 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
+                         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 /) )                                                                     
                                                                                                                                       
@@ -764,11 +767,11 @@
  
  
-    USE control_parameters,                                                 &
+    USE control_parameters,                                                    &
         ONLY:  data_output, message_string
 
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  unit     !< 
-    CHARACTER (LEN=*) ::  var !<
+    CHARACTER (LEN=*) ::  unit  !< 
+    CHARACTER (LEN=*) ::  var   !<
 
     INTEGER(iwp) :: i
@@ -780,5 +783,5 @@
        CASE ( 'm_soil' )
           IF (  .NOT.  land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                       'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
@@ -788,5 +791,5 @@
        CASE ( 't_soil' )
           IF (  .NOT.  land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                       'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
@@ -794,85 +797,85 @@
           unit = 'K'   
              
-       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf*', 'm_liq*',&
-              'qsws_liq*', 'qsws_soil*', 'qsws_veg*',  &
+       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf*', 'm_liq*',         &
+              'qsws_liq*', 'qsws_soil*', 'qsws_veg*',                          &
               'r_a*', 'r_s*' )
           IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
-             message_string = 'illegal value for data_output: "' //         &
-                              TRIM( var ) // '" & only 2d-horizontal ' //   &
+             message_string = 'illegal value for data_output: "' //            &
+                              TRIM( var ) // '" & only 2d-horizontal ' //      &
                               'cross sections are allowed for this value'
              CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
           ENDIF
           IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
           IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
           IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
           IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
           IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
           IF ( TRIM( var ) == 'm_liq*'  .AND.  .NOT.  land_surface )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'qsws_liq*'  .AND.  .NOT. land_surface )   &
+          IF ( TRIM( var ) == 'qsws_liq*'  .AND.  .NOT. land_surface )         &
           THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'qsws_soil*'  .AND.  .NOT.  land_surface ) &
+          IF ( TRIM( var ) == 'qsws_soil*'  .AND.  .NOT.  land_surface )       &
           THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'qsws_veg*'  .AND.  .NOT. land_surface )   &
+          IF ( TRIM( var ) == 'qsws_veg*'  .AND.  .NOT. land_surface )         &
           THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
+          IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface )             &
           THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
+          IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface )             &
           THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                               'res land_surface = .TRUE.'
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
 
-          IF ( TRIM( var ) == 'lai*'   )  unit = 'none' 
-          IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
-          IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
-          IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
-          IF ( TRIM( var ) == 'ghf*')  unit = 'W/m2'
+          IF ( TRIM( var ) == 'lai*'   )      unit = 'none' 
+          IF ( TRIM( var ) == 'c_liq*' )      unit = 'none'
+          IF ( TRIM( var ) == 'c_soil*')      unit = 'none'
+          IF ( TRIM( var ) == 'c_veg*' )      unit = 'none'
+          IF ( TRIM( var ) == 'ghf*')         unit = 'W/m2'
           IF ( TRIM( var ) == 'm_liq*'     )  unit = 'm'
-          IF ( TRIM( var ) == 'qsws_liq*'  ) unit = 'W/m2'
-          IF ( TRIM( var ) == 'qsws_soil*' ) unit = 'W/m2'
-          IF ( TRIM( var ) == 'qsws_veg*'  ) unit = 'W/m2'
-          IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
-          IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
+          IF ( TRIM( var ) == 'qsws_liq*'  )  unit = 'W/m2'
+          IF ( TRIM( var ) == 'qsws_soil*' )  unit = 'W/m2'
+          IF ( TRIM( var ) == 'qsws_veg*'  )  unit = 'W/m2'
+          IF ( TRIM( var ) == 'r_a*')         unit = 's/m'     
+          IF ( TRIM( var ) == 'r_s*')         unit = 's/m' 
              
        CASE DEFAULT
@@ -892,5 +895,5 @@
  SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
  
-    USE control_parameters,                                                 &
+    USE control_parameters,                                                    &
         ONLY:  data_output_pr, message_string
 
@@ -914,6 +917,6 @@
        CASE ( 't_soil', '#t_soil' )
           IF (  .NOT.  land_surface )  THEN
-             message_string = 'data_output_pr = ' //                        &
-                              TRIM( data_output_pr(var_count) ) // ' is' // &
+             message_string = 'data_output_pr = ' //                           &
+                              TRIM( data_output_pr(var_count) ) // ' is' //    &
                               'not implemented for land_surface = .FALSE.'
              CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
@@ -932,6 +935,6 @@
        CASE ( 'm_soil', '#m_soil' )
           IF (  .NOT.  land_surface )  THEN
-             message_string = 'data_output_pr = ' //                        &
-                              TRIM( data_output_pr(var_count) ) // ' is' // &
+             message_string = 'data_output_pr = ' //                           &
+                              TRIM( data_output_pr(var_count) ) // ' is' //    &
                               ' not implemented for land_surface = .FALSE.'
              CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
@@ -983,6 +986,6 @@
          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 )
@@ -1010,6 +1013,6 @@
        IF ( vegetation_type == 0 )  THEN
           IF ( min_canopy_resistance == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user defined)'//              &
-                              'requires setting of min_canopy_resistance'//       &
+             message_string = 'vegetation_type = 0 (user defined)'//           &
+                              'requires setting of min_canopy_resistance'//    &
                               '/= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1017,6 +1020,6 @@
 
           IF ( leaf_area_index == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of leaf_area_index'//             &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of leaf_area_index'//          &
                               '/= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1024,6 +1027,6 @@
 
           IF ( vegetation_coverage == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of vegetation_coverage'//         &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of vegetation_coverage'//      &
                               '/= 9999999.9'
                 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1031,6 +1034,6 @@
 
           IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of'//                             &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of'//                          &
                               'canopy_resistance_coefficient /= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1038,6 +1041,6 @@
 
           IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of lambda_surface_stable'//       &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of lambda_surface_stable'//    &
                               '/= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1045,6 +1048,6 @@
 
           IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of lambda_surface_unstable'//     &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of lambda_surface_unstable'//  &
                               '/= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1052,6 +1055,6 @@
 
           IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of f_shortwave_incoming'//        &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of f_shortwave_incoming'//     &
                               '/= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1059,6 +1062,6 @@
 
           IF ( z0_vegetation == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of z0_vegetation'//               &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of z0_vegetation'//            &
                               '/= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
@@ -1066,11 +1069,19 @@
 
           IF ( z0h_vegetation == 9999999.9_wp )  THEN
-             message_string = 'vegetation_type = 0 (user_defined)'//              &
-                              'requires setting of z0h_vegetation'//              &
+             message_string = 'vegetation_type = 0 (user_defined)'//           &
+                              'requires setting of z0h_vegetation'//           &
                               '/= 9999999.9'
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
           ENDIF
        ENDIF
-       
+     
+       IF ( vegetation_type == 1 )  THEN
+          IF ( vegetation_coverage /= 9999999.9_wp  .AND.  vegetation_coverage &
+               /= 0.0_wp )  THEN
+               message_string = 'vegetation_type = 1 (bare soil)'//            &
+                              ' requires vegetation_coverage = 0'//            &
+             CALL message( 'check_parameters', 'PA0471', 1, 2, 0, 6, 0 )
+          ENDIF
+  
     ENDIF
     
@@ -1165,6 +1176,6 @@
 
        IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of alpha_vangenuchten'//       &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of alpha_vangenuchten'//          &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1172,6 +1183,6 @@
 
        IF ( l_vangenuchten == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of l_vangenuchten'//           &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of l_vangenuchten'//              &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1179,6 +1190,6 @@
 
        IF ( n_vangenuchten == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of n_vangenuchten'//           &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of n_vangenuchten'//              &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1186,6 +1197,6 @@
 
        IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of hydraulic_conductivity'//   &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of hydraulic_conductivity'//      &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1193,6 +1204,6 @@
 
        IF ( saturation_moisture == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of saturation_moisture'//      &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of saturation_moisture'//         &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1200,6 +1211,6 @@
 
        IF ( field_capacity == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of field_capacity'//           &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of field_capacity'//              &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1207,6 +1218,6 @@
 
        IF ( wilting_point == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of wilting_point'//            &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of wilting_point'//               &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1214,6 +1225,6 @@
 
        IF ( residual_moisture == 9999999.9_wp )  THEN
-          message_string = 'soil_type = 0 (user_defined)'//                 &
-                           'requires setting of residual_moisture'//        &
+          message_string = 'soil_type = 0 (user_defined)'//                    &
+                           'requires setting of residual_moisture'//           &
                            '/= 9999999.9'
           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
@@ -1255,5 +1266,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 )
@@ -1286,5 +1297,5 @@
        dz_soil_center(k) = zs(k+1) - zs(k)
        IF ( dz_soil_center(k) == 0.0_wp )  THEN
-          message_string = 'invalid soil layer configuration found ' //     &
+          message_string = 'invalid soil layer configuration found ' //        &
                            '(dz_soil_center(k) = 0.0)'
           CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
@@ -1363,11 +1374,10 @@
     IF ( horizontal )  THEN
        surf              => surf_lsm_h
-
        surf_t_surface    => t_surface_h
        surf_t_surface_p  => t_surface_h_p
        surf_tt_surface_m => tt_surface_h_m
-       surf_m_liq     => m_liq_h
-       surf_m_liq_p   => m_liq_h_p
-       surf_tm_liq_m  => tm_liq_h_m
+       surf_m_liq        => m_liq_h
+       surf_m_liq_p      => m_liq_h_p
+       surf_tm_liq_m     => tm_liq_h_m
        surf_m_soil       => m_soil_h
        surf_t_soil       => t_soil_h
@@ -1378,11 +1388,10 @@
     ELSE
        surf              => surf_lsm_v(l)
-
        surf_t_surface    => t_surface_v(l)
        surf_t_surface_p  => t_surface_v_p(l)
        surf_tt_surface_m => tt_surface_v_m(l)
-       surf_m_liq     => m_liq_v(l)
-       surf_m_liq_p   => m_liq_v_p(l)
-       surf_tm_liq_m  => tm_liq_v_m(l)
+       surf_m_liq        => m_liq_v(l)
+       surf_m_liq_p      => m_liq_v_p(l)
+       surf_tm_liq_m     => tm_liq_v_m(l)
        surf_m_soil       => m_soil_v(l)
        surf_t_soil       => t_soil_v(l)
@@ -1519,7 +1528,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 * 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)) )
 !
 !--    f2: correction for soil moisture availability to plants (the 
@@ -1571,5 +1580,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     )
@@ -1586,5 +1596,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)              &
@@ -1706,6 +1716,6 @@
 !--    often reached, when no oscillations would occur (causes immense
 !--    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
+       IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )       &
+            > 0.2_wp  .AND. unscheduled_radiation_calls )  THEN
           force_radiation_call_l = .TRUE.
        ENDIF
@@ -1725,5 +1735,5 @@
                      ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
 
-       surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)          &
+       surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)             &
                                              - surf_t_soil%var_2d(nzb_soil,m) )
 
@@ -1732,17 +1742,17 @@
 
        IF ( humidity )  THEN
-          surf%qsws(m)  = - f_qsws * ( qv1 - q_s + dq_s_dt                  &
+          surf%qsws(m)  = - f_qsws * ( qv1 - 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  * ( qv1 - 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 * ( qv1 - 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  * ( qv1 - q_s                      &
                               + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
                               * surf_t_surface_p%var_1d(m) )
@@ -1776,10 +1786,10 @@
 !--          as runoff as qsws_soil is then not used in the soil model
              IF ( surf_m_liq%var_1d(m) /= m_liq_max )  THEN
-                surf%qsws_liq(m) = surf%qsws_liq(m)                      &
+                surf%qsws_liq(m) = surf%qsws_liq(m)                            &
                                  + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
                                  * hyrho(k+k_off)                              &
                                  * 0.001_wp * rho_l * l_v
              ELSE
-                surf%qsws_soil(m) = surf%qsws_soil(m)                    &
+                surf%qsws_soil(m) = surf%qsws_soil(m)                          &
                                  + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
                                  * hyrho(k+k_off)                              &
@@ -1789,5 +1799,5 @@
 !--          Add precipitation to bare soil according to the bare soil
 !--          coverage.
-             surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp            &
+             surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp                  &
                                - surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
                                * hyrho(k+k_off)                                &
@@ -1822,5 +1832,5 @@
  
           tend = - surf%qsws_liq(m) * drho_l_lv
-          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *        &
+          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *              &
                                         ( tsc(2) * tend +                      &
                                           tsc(3) * surf_tm_liq_m%var_1d(m) )
@@ -1841,5 +1851,5 @@
              ELSEIF ( intermediate_timestep_count <                            &
                       intermediate_timestep_count_max )  THEN
-                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +               &
+                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +                  &
                                               5.3125_wp * surf_tm_liq_m%var_1d(m)
              ENDIF
@@ -2404,9 +2414,9 @@
                 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)    &
+                                                 * 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)     &
+                                                 * dz_soil(nzb_soil)           &
                                                  * 0.25_wp
                 surf_lsm_h%lambda_h_def(m)    = pavement_heat_conduct                                     
@@ -2455,6 +2465,8 @@
 
              surf_lsm_h%us(m)   = 0.1_wp
-             surf_lsm_h%ts(m)   = ( pt1 - 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
+             surf_lsm_h%ts(m)   = ( pt1 - 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
@@ -2491,10 +2503,10 @@
                    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) &
+                   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)      &
+                                                * dz_soil(nzb_soil)            &
                                                 * 0.25_wp
                    surf_lsm_v(l)%lambda_h_def(m)   = pavement_heat_conduct                                     
@@ -2765,6 +2777,6 @@
 !
 !--    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 ( 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)                  )
@@ -2777,6 +2789,6 @@
 !--    Vertical surfaces
        DO  l = 0, 3
-          ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns)                     )
-          ALLOCATE ( m_liq_v_p(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
+          ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns)                        )
+          ALLOCATE ( m_liq_v_p(l)%var_1d(1:surf_lsm_v(l)%ns)                      )
           ALLOCATE ( t_surface_v(l)%var_1d(1:surf_lsm_v(l)%ns)                    )
           ALLOCATE ( t_surface_v_p(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
@@ -2789,6 +2801,6 @@
 !
 !--    Horizontal surfaces
-       ALLOCATE ( m_liq_h_1%var_1d(1:surf_lsm_h%ns)                   )
-       ALLOCATE ( m_liq_h_2%var_1d(1:surf_lsm_h%ns)                   )
+       ALLOCATE ( m_liq_h_1%var_1d(1:surf_lsm_h%ns)                      )
+       ALLOCATE ( m_liq_h_2%var_1d(1:surf_lsm_h%ns)                      )
        ALLOCATE ( t_surface_h_1%var_1d(1:surf_lsm_h%ns)                  )
        ALLOCATE ( t_surface_h_2%var_1d(1:surf_lsm_h%ns)                  )
@@ -2800,6 +2812,6 @@
 !--    Vertical surfaces
        DO  l = 0, 3
-          ALLOCATE ( m_liq_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
-          ALLOCATE ( m_liq_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
+          ALLOCATE ( m_liq_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                      )
+          ALLOCATE ( m_liq_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                      )
           ALLOCATE ( t_surface_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
           ALLOCATE ( t_surface_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
@@ -2814,5 +2826,5 @@
 !--    Allocate intermediate timestep arrays
 !--    Horizontal surfaces
-       ALLOCATE ( tm_liq_h_m%var_1d(1:surf_lsm_h%ns)                  )
+       ALLOCATE ( tm_liq_h_m%var_1d(1:surf_lsm_h%ns)                     )
        ALLOCATE ( tt_surface_h_m%var_1d(1:surf_lsm_h%ns)                 )
        ALLOCATE ( tm_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
@@ -2821,5 +2833,5 @@
 !--    Horizontal surfaces
        DO  l = 0, 3
-          ALLOCATE ( tm_liq_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
+          ALLOCATE ( tm_liq_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                     )
           ALLOCATE ( tt_surface_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                 )
           ALLOCATE ( tm_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
@@ -2873,5 +2885,5 @@
           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)%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)              )
@@ -2907,5 +2919,5 @@
        t_surface_h => t_surface_h_1; t_surface_h_p => t_surface_h_2
        m_soil_h    => m_soil_h_1;    m_soil_h_p    => m_soil_h_2
-       m_liq_h  => m_liq_h_1;  m_liq_h_p  => m_liq_h_2
+       m_liq_h     => m_liq_h_1;     m_liq_h_p     => m_liq_h_2
 !
 !--    Vertical surfaces
@@ -2913,5 +2925,5 @@
        t_surface_v => t_surface_v_1; t_surface_v_p => t_surface_v_2
        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
+       m_liq_v     => m_liq_v_1;     m_liq_v_p     => m_liq_v_2
 #endif
 
@@ -3029,5 +3041,4 @@
        IF ( horizontal )  THEN
           surf           => surf_lsm_h
-
           surf_m_soil    => m_soil_h
           surf_m_soil_p  => m_soil_h_p
@@ -3038,5 +3049,4 @@
        ELSE
           surf           => surf_lsm_v(l)
-
           surf_m_soil    => m_soil_v(l)
           surf_m_soil_p  => m_soil_v_p(l)
@@ -3061,6 +3071,6 @@
 !--                Calculate volumetric heat capacity of the soil, taking 
 !--                into account water content 
-                   surf%rho_c_total(k,m) = (rho_c_soil *                 &
-                                               ( 1.0_wp - surf%m_sat(k,m) )&
+                   surf%rho_c_total(k,m) = (rho_c_soil *                       &
+                                               ( 1.0_wp - surf%m_sat(k,m) )    &
                                                + rho_c_water * surf_m_soil%var_2d(k,m) )
 
@@ -3068,5 +3078,5 @@
 !--                Calculate soil heat conductivity at the center of the soil
 !--                layers
-                   lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(k,m)) *      &
+                   lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(k,m)) *    &
                                   lambda_h_water ** surf_m_soil%var_2d(k,m)
 
@@ -3087,8 +3097,8 @@
                      .AND.  zs(k+1) > pavement_depth )                         &
                 THEN
-                   surf%lambda_h(k,m) = ( pavement_depth - zs(k) ) * ddz_soil_center(k+1)   &
-                                      * lambda_temp(k)                               &
-                                      + ( zs(k+1) - pavement_depth ) * ddz_soil_center(k+1) &
-                                      * lambda_temp(k+1)
+                   surf%lambda_h(k,m) = ( pavement_depth - zs(k) )             &
+                                      * ddz_soil_center(k+1) * lambda_temp(k)  &
+                                      + ( zs(k+1) - pavement_depth )           &
+                                      * ddz_soil_center(k+1) * lambda_temp(k+1)
                 ELSE
                    surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
@@ -3103,26 +3113,26 @@
              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)
-
-             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) )
+                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) )
 
 !
@@ -3136,6 +3146,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
@@ -3148,30 +3159,27 @@
 !--             Calculate soil diffusivity at the center of the soil layers
                 lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat     &
-                                  / surf%m_sat(k,m) ) * ( MAX( surf_m_soil%var_2d(k,m),    &
-                                  surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**(&
-                                                            b_ch + 2.0_wp )
+                                 / surf%m_sat(k,m) ) * (                       &
+                                 MAX( surf_m_soil%var_2d(k,m),                 &
+                                 surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**(     &
+                                 b_ch + 2.0_wp )
 
 !
 !--             Parametrization of Van Genuchten
 !--             Calculate the hydraulic conductivity after Van Genuchten (1980)
-                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) )       &
-                           )                                                &
-                         )**(                                               &
-                       surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) &
-                            ) - 1.0_wp                                      &
+                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) )  &
+                           )                                                   &
+                         )**( 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)
 
-                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) ) )
+                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) ) )
 
              ENDDO
@@ -3198,7 +3206,7 @@
                 DO  k = nzb_soil, nzt_soil-1
                     
-                   surf%lambda_w(k,m) = ( lambda_temp(k+1) +             &
+                   surf%lambda_w(k,m) = ( lambda_temp(k+1) +                   &
                                                 lambda_temp(k) ) * 0.5_wp
-                   surf%gamma_w(k,m)  = ( gamma_temp(k+1)  +             &
+                   surf%gamma_w(k,m)  = ( gamma_temp(k+1)  +                   &
                                                 gamma_temp(k) )  * 0.5_wp
 
@@ -3232,5 +3240,6 @@
                 DO  k = nzb_soil, nzt_soil
                     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)
+                       m_total = m_total + surf%root_fr(k,m)                   &
+                                 * surf_m_soil%var_2d(k,m)
                     ENDIF
                 ENDDO  
@@ -3238,6 +3247,6 @@
                    DO  k = nzb_soil, nzt_soil
                       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
+                         root_extr(k) = surf%root_fr(k,m)                      &
+                                        * surf_m_soil%var_2d(k,m) / m_total
                       ELSE
                          root_extr(k) = 0.0_wp
@@ -3249,35 +3258,36 @@
                 tend(:) = 0.0_wp
 
-                tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (       &
-                         surf_m_soil%var_2d(nzb_soil+1,m) - surf_m_soil%var_2d(nzb_soil,m) )           &
-                         * ddz_soil_center(nzb_soil) - surf%gamma_w(nzb_soil,m) - &
-                                                     (                         &
-                            root_extr(nzb_soil) * surf%qsws_veg(m)    &
-                            + surf%qsws_soil(m) ) * drho_l_lv )       &
+                tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (             &
+                         surf_m_soil%var_2d(nzb_soil+1,m)                      &
+                         - surf_m_soil%var_2d(nzb_soil,m) )                    &
+                         * ddz_soil_center(nzb_soil) - surf%gamma_w(nzb_soil,m)&
+                         - ( root_extr(nzb_soil) * surf%qsws_veg(m)            &
+                            + surf%qsws_soil(m) ) * drho_l_lv )                &
                             * ddz_soil(nzb_soil)
 
                 DO  k = nzb_soil+1, nzt_soil-1
-                   tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)      &
-                             - surf_m_soil%var_2d(k,m) ) * ddz_soil_center(k)                   &
-                             - surf%gamma_w(k,m)                         &
-                             - surf%lambda_w(k-1,m) * ( surf_m_soil%var_2d(k,m) -    &
-                             surf_m_soil%var_2d(k-1,m)) * ddz_soil_center(k-1)                      &
-                             + surf%gamma_w(k-1,m) - (root_extr(k)       &
-                             * surf%qsws_veg(m) * drho_l_lv)          &
+                   tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)  &
+                             - surf_m_soil%var_2d(k,m) ) * ddz_soil_center(k)    &
+                             - surf%gamma_w(k,m)                                 &
+                             - surf%lambda_w(k-1,m) * ( surf_m_soil%var_2d(k,m)  &
+                             - surf_m_soil%var_2d(k-1,m)) * ddz_soil_center(k-1) &
+                             + surf%gamma_w(k-1,m) - (root_extr(k)               &
+                             * surf%qsws_veg(m) * drho_l_lv)                     &
                              ) * ddz_soil(k)
                 ENDDO
-                tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)            &
-                                     - surf%lambda_w(nzt_soil-1,m)       &
-                                     * ( surf_m_soil%var_2d(nzt_soil,m)                    &
-                                     - surf_m_soil%var_2d(nzt_soil-1,m))                   &
-                                     * ddz_soil_center(nzt_soil-1)                      &
-                                     + surf%gamma_w(nzt_soil-1,m) - (    &
-                                       root_extr(nzt_soil)                     &
-                                     * surf%qsws_veg(m) * drho_l_lv ) &
+                tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)                  &
+                                   - surf%lambda_w(nzt_soil-1,m)               &
+                                   * ( surf_m_soil%var_2d(nzt_soil,m)          &
+                                   - surf_m_soil%var_2d(nzt_soil-1,m))         &
+                                   * ddz_soil_center(nzt_soil-1)               &
+                                   + surf%gamma_w(nzt_soil-1,m) - (            &
+                                   root_extr(nzt_soil)                         &
+                                   * surf%qsws_veg(m) * drho_l_lv )            &
                                   ) * ddz_soil(nzt_soil)             
 
-                surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_m_soil%var_2d(nzb_soil:nzt_soil,m)    &
-                                                + dt_3d * ( tsc(2) * tend(:)   &
-                                                + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
+                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) )   
    
 !
@@ -3297,6 +3307,7 @@
                             intermediate_timestep_count_max )  THEN
                       DO  k = nzb_soil, nzt_soil
-                         surf_tm_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) + 5.3125_wp     &
-                                  * surf_tm_soil_m%var_2d(k,m)
+                         surf_tm_soil_m%var_2d(k,m) = -9.5625_wp * tend(k)     &
+                                                    + 5.3125_wp                &
+                                                    * surf_tm_soil_m%var_2d(k,m)
                       ENDDO
                    ENDIF
@@ -3350,6 +3361,6 @@
              t_soil_h     => t_soil_h_1;    t_soil_h_p     => t_soil_h_2
              IF ( humidity )  THEN
-                m_soil_h    => m_soil_h_1;   m_soil_h_p    => m_soil_h_2
-                m_liq_h  => m_liq_h_1; m_liq_h_p  => m_liq_h_2
+                m_soil_h  => m_soil_h_1;    m_soil_h_p     => m_soil_h_2
+                m_liq_h   => m_liq_h_1;     m_liq_h_p      => m_liq_h_2
              ENDIF
 !
@@ -3358,6 +3369,6 @@
              t_soil_v     => t_soil_v_1;    t_soil_v_p     => t_soil_v_2
              IF ( humidity )  THEN
-                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
+                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
 
@@ -3370,6 +3381,6 @@
              t_soil_h     => t_soil_h_2;    t_soil_h_p     => t_soil_h_1
              IF ( humidity )  THEN
-                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
+                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
 !
@@ -3378,6 +3389,6 @@
              t_soil_v     => t_soil_v_2;    t_soil_v_p     => t_soil_v_1
              IF ( humidity )  THEN
-                m_soil_v    => m_soil_v_2;   m_soil_v_p    => m_soil_v_1
-                m_liq_v  => m_liq_v_2; m_liq_v_p  => m_liq_v_1
+                m_soil_v  => m_soil_v_2;    m_soil_v_p     => m_soil_v_1
+                m_liq_v   => m_liq_v_2;     m_liq_v_p      => m_liq_v_1
              ENDIF
 
@@ -3561,5 +3572,5 @@
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                qsws_liq_av(j,i) = qsws_liq_av(j,i) +                    &
+                qsws_liq_av(j,i) = qsws_liq_av(j,i) +                          &
                                       surf_lsm_h%qsws_liq(m)
              ENDDO
@@ -3569,5 +3580,5 @@
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                qsws_soil_av(j,i) = qsws_soil_av(j,i) +                  &
+                qsws_soil_av(j,i) = qsws_soil_av(j,i) +                        &
                                        surf_lsm_h%qsws_soil(m)
              ENDDO
@@ -3577,5 +3588,5 @@
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                qsws_veg_av(j,i) = qsws_veg_av(j,i) +                    &
+                qsws_veg_av(j,i) = qsws_veg_av(j,i) +                          &
                                       surf_lsm_h%qsws_veg(m)
              ENDDO
@@ -3616,5 +3627,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   c_liq_av(j,i) = c_liq_av(j,i)                               &
+                                   / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3623,5 +3635,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   c_soil_av(j,i) = c_soil_av(j,i)                             &
+                                    / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3630,5 +3643,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   c_veg_av(j,i) = c_veg_av(j,i)                               &
+                                   / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3637,5 +3651,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   ghf_av(j,i) = ghf_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   ghf_av(j,i) = ghf_av(j,i)                                   &
+                                 / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3644,5 +3659,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   lai_av(j,i) = lai_av(j,i)                                   &
+                                 / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3651,5 +3667,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   m_liq_av(j,i) = m_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   m_liq_av(j,i) = m_liq_av(j,i)                               &
+                                   / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3659,5 +3676,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_soil, nzt_soil
-                      m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      m_soil_av(k,j,i) = m_soil_av(k,j,i)                      &
+                                         / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -3667,5 +3685,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   qsws_liq_av(j,i) = qsws_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   qsws_liq_av(j,i) = qsws_liq_av(j,i)                         &
+                                      / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3674,5 +3693,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   qsws_soil_av(j,i) = qsws_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   qsws_soil_av(j,i) = qsws_soil_av(j,i)                       &
+                                       / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3681,5 +3701,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   qsws_veg_av(j,i) = qsws_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   qsws_veg_av(j,i) = qsws_veg_av(j,i)                         &
+                                      / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -3703,5 +3724,6 @@
                 DO  j = nys, nyn
                    DO  k = nzb_soil, nzt_soil
-                      t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      t_soil_av(k,j,i) = t_soil_av(k,j,i)                      &
+                                         / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -4239,25 +4261,25 @@
 
     REAL(wp),                                                                  &
-       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: &
           tmp_2d   !< 
 
     REAL(wp),                                                                  &
-       DIMENSION(nzb_soil:nzt_soil+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+       DIMENSION(nzb_soil:nzt_soil+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: &
           tmp_3d   !< 
 
     REAL(wp),                                                                  &
-       DIMENSION(nzb_soil:nzt_soil,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+       DIMENSION(nzb_soil:nzt_soil,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: &
           tmp_3d2   !< 
 
     REAL(wp),                                                                  &
-       DIMENSION(1:surf_lsm_h%ns) ::                                  &
+       DIMENSION(1:surf_lsm_h%ns) ::                                           &
           tmp_walltype_1d   !< 
 
     REAL(wp),                                                                  &
-       DIMENSION(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ::              &
+       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) ::                &
+       DIMENSION(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ::                         &
           tmp_walltype_2d2  !< 
 
@@ -4332,5 +4354,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   m_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
+                   m_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -4354,5 +4376,5 @@
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
+                   qsws_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =     &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
                 CASE ( 'qsws_soil_av' )
@@ -4361,5 +4383,5 @@
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                   qsws_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =    &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -4369,5 +4391,5 @@
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
+                   qsws_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =     &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
