Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 2269)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 2270)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Revised numbering (removed 2 timeseries)
+! 
+! 2259 2017-06-08 09:09:11Z gronemeier
 ! Implemented synthetic turbulence generator
 !
@@ -923,4 +926,11 @@
     ENDIF
 
+
+!
+!-- When the land surface model is used, the flux output must be dynamic.
+    IF ( land_surface )  THEN
+       flux_output_mode = 'dynamic'
+    ENDIF
+    
 !
 !-- set the flux output units according to flux_output_mode
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 2269)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 2270)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Revised numbering (removed 2 timeseries)
+! 
+! 2252 2017-06-07 09:35:37Z knoop
 ! perturbation pressure now depending on flux_output_mode
 ! 
@@ -399,5 +402,5 @@
              IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
                                            * waterflux_output_conversion  ! w*q*
-             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
+             IF ( passive_scalar ) sums_l(:,114,i) = sums_wsss_ws_l(:,i)  ! w*s*
           ENDDO
 
@@ -484,5 +487,5 @@
              DO  j =  nys, nyn
                 DO  k = nzb, nzt+1
-                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i)             &
+                   sums_l(k,115,tn)  = sums_l(k,115,tn) + s(k,j,i)             &
                                     * rmask(j,i,sr)                            &
                                     * MERGE( 1.0_wp, 0.0_wp,                   &
@@ -512,5 +515,5 @@
              ENDIF
              IF ( passive_scalar )  THEN
-                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
+                sums_l(:,115,0) = sums_l(:,115,0) + sums_l(:,115,i)
              ENDIF
           ENDDO
@@ -553,5 +556,5 @@
        IF ( passive_scalar )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,       &
+          CALL MPI_ALLREDUCE( sums_l(nzb,115,0), sums(nzb,115), nzt+2-nzb,       &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
        ENDIF
@@ -569,5 +572,5 @@
           ENDIF
        ENDIF
-       IF ( passive_scalar )  sums(:,117) = sums_l(:,117,0)
+       IF ( passive_scalar )  sums(:,115) = sums_l(:,115,0)
 #endif
 
@@ -607,5 +610,5 @@
 !
 !--    Passive scalar
-       IF ( passive_scalar )  hom(:,1,117,sr) = sums(:,117) /                  &
+       IF ( passive_scalar )  hom(:,1,115,sr) = sums(:,115) /                  &
             ngp_2dh_s_inner(:,sr)                    ! s 
 
@@ -654,6 +657,6 @@
                 ENDIF
                 IF ( passive_scalar )  THEN
-                   sums_l(k,118,tn) = sums_l(k,118,tn) + &
-                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)&
+                   sums_l(k,116,tn) = sums_l(k,116,tn) + &
+                                  ( s(k,j,i)-hom(k,1,115,sr) )**2 * rmask(j,i,sr)&
                                                                   * flag
                 ENDIF
@@ -889,5 +892,5 @@
 !--             Passive scalar flux
                 IF ( passive_scalar )  THEN
-                   sums_l(k,119,tn) = sums_l(k,119,tn)                         &
+                   sums_l(k,117,tn) = sums_l(k,117,tn)                         &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                   * ( s(k+1,j,i) - s(k,j,i) )  &
@@ -949,5 +952,5 @@
                          ENDIF
                          IF ( passive_scalar )  THEN
-                            sums_l(k+ki,119,tn) = sums_l(k+ki,119,tn) +                     &
+                            sums_l(k+ki,117,tn) = sums_l(k+ki,117,tn) +                     &
                                         surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
                          ENDIF
@@ -1001,5 +1004,5 @@
                    ENDIF
                    IF ( passive_scalar )  THEN
-                      sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
+                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
                                         surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
                    ENDIF
@@ -1051,5 +1054,5 @@
                    ENDIF
                    IF ( passive_scalar )  THEN
-                      sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
+                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
                                         surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
                    ENDIF
@@ -1063,15 +1066,15 @@
                 IF ( surf_def_h(0)%ns >= 1 )  THEN
                    m = surf_def_h(0)%start_index(j,i)
-                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
+                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                      &
                                         surf_def_h(0)%ol(m)  * rmask(j,i,sr) ! L
                 ENDIF
                 IF ( surf_lsm_h%ns >= 1 )  THEN
                    m = surf_lsm_h%start_index(j,i)
-                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
+                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                      &
                                         surf_lsm_h%ol(m)  * rmask(j,i,sr) ! L
                 ENDIF
                 IF ( surf_usm_h%ns >= 1 )  THEN
                    m = surf_usm_h%start_index(j,i)
-                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
+                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                      &
                                         surf_usm_h%ol(m)  * rmask(j,i,sr) ! L
                 ENDIF
@@ -1079,16 +1082,16 @@
 
              IF ( radiation  .AND.  radiation_scheme /= 'constant' )  THEN
-                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_net(j,i)
-                sums_l(nzb,102,tn)  = sums_l(nzb,102,tn)  + rad_lw_in(nzb,j,i)
-                sums_l(nzb,103,tn)  = sums_l(nzb,103,tn)  + rad_lw_out(nzb,j,i)
-                sums_l(nzb,104,tn)  = sums_l(nzb,104,tn)  + rad_sw_in(nzb,j,i)
-                sums_l(nzb,105,tn)  = sums_l(nzb,105,tn)  + rad_sw_out(nzb,j,i)
+                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  + rad_net(j,i)
+                sums_l(nzb,100,tn)  = sums_l(nzb,100,tn)  + rad_lw_in(nzb,j,i)
+                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_lw_out(nzb,j,i)
+                sums_l(nzb,102,tn)  = sums_l(nzb,102,tn)  + rad_sw_in(nzb,j,i)
+                sums_l(nzb,103,tn)  = sums_l(nzb,103,tn)  + rad_sw_out(nzb,j,i)
 
 #if defined ( __rrtmg )
                 IF ( radiation_scheme == 'rrtmg' )  THEN
-                   sums_l(nzb,110,tn)  = sums_l(nzb,110,tn)  + rrtm_aldif(0,j,i)
-                   sums_l(nzb,111,tn)  = sums_l(nzb,111,tn)  + rrtm_aldir(0,j,i)
-                   sums_l(nzb,112,tn)  = sums_l(nzb,112,tn)  + rrtm_asdif(0,j,i)
-                   sums_l(nzb,113,tn)  = sums_l(nzb,113,tn)  + rrtm_asdir(0,j,i)
+                   sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  + rrtm_aldif(0,j,i)
+                   sums_l(nzb,109,tn)  = sums_l(nzb,109,tn)  + rrtm_aldir(0,j,i)
+                   sums_l(nzb,110,tn)  = sums_l(nzb,110,tn)  + rrtm_asdif(0,j,i)
+                   sums_l(nzb,111,tn)  = sums_l(nzb,111,tn)  + rrtm_asdir(0,j,i)
                 ENDIF
 #endif
@@ -1144,5 +1147,5 @@
                 ENDIF
                 IF ( passive_scalar )  THEN
-                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
+                   sums_l(nzt,117,tn) = sums_l(nzt,117,tn) + &
                                         surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
                 ENDIF
@@ -1247,7 +1250,7 @@
                 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
                      .OR. sr /= 0 ) )  THEN
-                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +             &
-                                    s(k+1,j,i) - hom(k+1,1,117,sr) )
-                   sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
+                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +             &
+                                    s(k+1,j,i) - hom(k+1,1,115,sr) )
+                   sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *      &
                                                        rmask(j,i,sr) * flag
                 ENDIF
@@ -1277,12 +1280,10 @@
              IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
                   j >= nys  .AND.  j <= nyn )  THEN 
-                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf_eb(m)
-                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%shf_eb(m)
-                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_eb(m)
-                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_liq_eb(m)
-                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%qsws_soil_eb(m)
-                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + surf_lsm_h%qsws_veg_eb(m)
-                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + surf_lsm_h%r_a(m)
-                sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ surf_lsm_h%r_s(m)
+                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)
+                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)
+                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_soil(m)
+                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_veg(m)
+                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%r_a(m)
+                sums_l(nzb,98,tn) = sums_l(nzb,98,tn)+ surf_lsm_h%r_s(m)
              ENDIF
           ENDDO
@@ -1379,7 +1380,7 @@
                    ENDIF
                    IF ( passive_scalar )  THEN
-                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +           &
-                                      s(k+1,j,i) - hom(k+1,1,117,sr) )
-                      sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *    &
+                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +           &
+                                      s(k+1,j,i) - hom(k+1,1,115,sr) )
+                      sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *    &
                                         rmask(j,i,sr) * flag
                    ENDIF
@@ -1581,19 +1582,19 @@
                    flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
-                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
+                   sums_l(k,100,tn)  = sums_l(k,100,tn)  + rad_lw_in(k,j,i)    &
                                        * rmask(j,i,sr) * flag
-                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
+                   sums_l(k,101,tn)  = sums_l(k,101,tn)  + rad_lw_out(k,j,i)   &
                                        * rmask(j,i,sr) * flag
-                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
+                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_sw_in(k,j,i)    &
                                        * rmask(j,i,sr) * flag
-                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
+                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_sw_out(k,j,i)   &
                                        * rmask(j,i,sr) * flag
-                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
+                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_lw_cs_hr(k,j,i) &
                                        * rmask(j,i,sr) * flag
-                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
+                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_lw_hr(k,j,i)    &
                                        * rmask(j,i,sr) * flag
-                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
+                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_sw_cs_hr(k,j,i) &
                                        * rmask(j,i,sr) * flag
-                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
+                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_sw_hr(k,j,i)    &
                                        * rmask(j,i,sr) * flag
                 ENDDO
@@ -1655,7 +1656,7 @@
           sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
           sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
-          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
-          sums(k,116)           = sums(k,116)           / ngp_2dh(sr)
-          sums(k,119)           = sums(k,119)           / ngp_2dh(sr)
+          sums(k,89:112)        = sums(k,89:112)        / ngp_2dh(sr)
+          sums(k,114)           = sums(k,114)           / ngp_2dh(sr)
+          sums(k,117)           = sums(k,117)           / ngp_2dh(sr)
           IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
              sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
@@ -1666,6 +1667,6 @@
              sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
              sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
-             sums(k,118)           = sums(k,118)           / ngp_2dh_s_inner(k,sr)
-             sums(k,120:pr_palm-2) = sums(k,120:pr_palm-2) / ngp_2dh_s_inner(k,sr)
+             sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
+             sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
           ENDIF
        ENDDO
@@ -1788,46 +1789,44 @@
           hom(:,1,91,sr) = sums(:,91)              ! m_soil
                                                    ! 92 is initial m_soil profile
-          hom(:,1,93,sr)  = sums(:,93)             ! ghf_eb
-          hom(:,1,94,sr)  = sums(:,94)             ! shf_eb
-          hom(:,1,95,sr)  = sums(:,95)             ! qsws_eb
-          hom(:,1,96,sr)  = sums(:,96)             ! qsws_liq_eb
-          hom(:,1,97,sr)  = sums(:,97)             ! qsws_soil_eb
-          hom(:,1,98,sr)  = sums(:,98)             ! qsws_veg_eb
-          hom(:,1,99,sr)  = sums(:,99)             ! r_a
-          hom(:,1,100,sr) = sums(:,100)            ! r_s
+          hom(:,1,93,sr)  = sums(:,93)             ! ghf
+          hom(:,1,94,sr)  = sums(:,94)             ! qsws_liq
+          hom(:,1,95,sr)  = sums(:,95)             ! qsws_soil
+          hom(:,1,96,sr)  = sums(:,96)             ! qsws_veg
+          hom(:,1,97,sr)  = sums(:,97)             ! r_a
+          hom(:,1,98,sr) = sums(:,98)              ! r_s
 
        ENDIF
 
        IF ( radiation )  THEN
-          hom(:,1,101,sr) = sums(:,101)            ! rad_net
-          hom(:,1,102,sr) = sums(:,102)            ! rad_lw_in
-          hom(:,1,103,sr) = sums(:,103)            ! rad_lw_out
-          hom(:,1,104,sr) = sums(:,104)            ! rad_sw_in
-          hom(:,1,105,sr) = sums(:,105)            ! rad_sw_out
+          hom(:,1,99,sr) = sums(:,99)            ! rad_net
+          hom(:,1,100,sr) = sums(:,100)            ! rad_lw_in
+          hom(:,1,101,sr) = sums(:,101)            ! rad_lw_out
+          hom(:,1,102,sr) = sums(:,102)            ! rad_sw_in
+          hom(:,1,103,sr) = sums(:,103)            ! rad_sw_out
 
           IF ( radiation_scheme == 'rrtmg' )  THEN
-             hom(:,1,106,sr) = sums(:,106)            ! rad_lw_cs_hr
-             hom(:,1,107,sr) = sums(:,107)            ! rad_lw_hr
-             hom(:,1,108,sr) = sums(:,108)            ! rad_sw_cs_hr
-             hom(:,1,109,sr) = sums(:,109)            ! rad_sw_hr
-
-             hom(:,1,110,sr) = sums(:,110)            ! rrtm_aldif
-             hom(:,1,111,sr) = sums(:,111)            ! rrtm_aldir
-             hom(:,1,112,sr) = sums(:,112)            ! rrtm_asdif
-             hom(:,1,113,sr) = sums(:,113)            ! rrtm_asdir
+             hom(:,1,104,sr) = sums(:,104)            ! rad_lw_cs_hr
+             hom(:,1,105,sr) = sums(:,105)            ! rad_lw_hr
+             hom(:,1,106,sr) = sums(:,106)            ! rad_sw_cs_hr
+             hom(:,1,107,sr) = sums(:,107)            ! rad_sw_hr
+
+             hom(:,1,108,sr) = sums(:,108)            ! rrtm_aldif
+             hom(:,1,109,sr) = sums(:,109)            ! rrtm_aldir
+             hom(:,1,110,sr) = sums(:,110)            ! rrtm_asdif
+             hom(:,1,111,sr) = sums(:,111)            ! rrtm_asdir
           ENDIF
        ENDIF
 
-       hom(:,1,114,sr) = sums(:,114)            !: L
+       hom(:,1,112,sr) = sums(:,112)            !: L
 
        IF ( passive_scalar )  THEN
-          hom(:,1,119,sr) = sums(:,119)     ! w"s"
-          hom(:,1,116,sr) = sums(:,116)     ! w*s*
-          hom(:,1,120,sr) = sums(:,119) + sums(:,116)    ! ws 
-          hom(:,1,118,sr) = sums(:,118)     ! s*2
-       ENDIF
-
-       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
-       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
+          hom(:,1,117,sr) = sums(:,117)     ! w"s"
+          hom(:,1,114,sr) = sums(:,114)     ! w*s*
+          hom(:,1,118,sr) = sums(:,117) + sums(:,114)    ! ws 
+          hom(:,1,116,sr) = sums(:,116)     ! s*2
+       ENDIF
+
+       hom(:,1,119,sr) = rho_air       ! rho_air in Kg/m^3
+       hom(:,1,120,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
 
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
@@ -1972,5 +1971,5 @@
 
        IF ( .NOT. neutral )  THEN
-          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
+          ts_value(22,sr) = hom(nzb,1,112,sr)          ! L
        ELSE
           ts_value(22,sr) = 1.0E10_wp
@@ -1980,5 +1979,5 @@
 
        IF ( passive_scalar )  THEN
-          ts_value(24,sr) = hom(nzb+13,1,119,sr)       ! w"s" ( to do ! )
+          ts_value(24,sr) = hom(nzb+13,1,117,sr)       ! w"s" ( to do ! )
           ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
        ENDIF
@@ -1987,27 +1986,25 @@
 !--    Collect land surface model timeseries
        IF ( land_surface )  THEN
-          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
-          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
-          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
-          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
-          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
-          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
-          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
-          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
+          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf
+          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! qsws_liq
+          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_soil
+          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_veg
+          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! r_a
+          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! r_s
        ENDIF
 !
 !--    Collect radiation model timeseries
        IF ( radiation )  THEN
-          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
-          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
-          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
-          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_sw_in
-          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
+          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
+          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_lw_in
+          ts_value(dots_rad+2,sr) = hom(nzb,1,101,sr)          ! rad_lw_out
+          ts_value(dots_rad+3,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
+          ts_value(dots_rad+4,sr) = hom(nzb,1,103,sr)          ! rad_sw_out
 
           IF ( radiation_scheme == 'rrtmg' )  THEN
-             ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr)          ! rrtm_aldif
-             ts_value(dots_rad+6,sr) = hom(nzb,1,111,sr)          ! rrtm_aldir
-             ts_value(dots_rad+7,sr) = hom(nzb,1,112,sr)          ! rrtm_asdif
-             ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr)          ! rrtm_asdir
+             ts_value(dots_rad+5,sr) = hom(nzb,1,108,sr)          ! rrtm_aldif
+             ts_value(dots_rad+6,sr) = hom(nzb,1,109,sr)          ! rrtm_aldir
+             ts_value(dots_rad+7,sr) = hom(nzb,1,110,sr)          ! rrtm_asdif
+             ts_value(dots_rad+8,sr) = hom(nzb,1,111,sr)          ! rrtm_asdir
           ENDIF
 
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 2269)
+++ /palm/trunk/SOURCE/header.f90	(revision 2270)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed Prandtl layer to constant flux layer
+! 
+! 2259 2017-06-08 09:09:11Z gronemeier
 ! Implemented synthetic turbulence generator
 !
@@ -2093,7 +2096,7 @@
 303 FORMAT (/' Bottom surface fluxes are used in diffusion terms at k=1')
 304 FORMAT (/' Top surface fluxes are used in diffusion terms at k=nzt')
-305 FORMAT (//'    Prandtl-Layer between bottom surface and first ', &
-               'computational u,v-level:'// &
-             '       zp = ',F6.2,' m   z0 =',F7.4,' m   z0h =',F8.5,&
+305 FORMAT (//'    Constant flux layer between bottom surface and first ',     &
+              'computational u,v-level:'// &
+             '       z_mo = ',F6.2,' m   z0 =',F7.4,' m   z0h =',F8.5,&
              ' m   kappa =',F5.2/ &
              '       Rif value range:   ',F8.2,' <= rif <=',F6.2)
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 2269)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 2270)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dots_num must be increased when LSM and/or radiation is used
+! 
+! 2259 2017-06-08 09:09:11Z gronemeier
 ! Implemented synthetic turbulence generator
 !
@@ -381,5 +384,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  radiation_init, radiation
+        ONLY:  radiation_init, radiation, radiation_scheme
     
     USE random_function_mod 
@@ -1863,4 +1866,16 @@
     ENDIF
 
+    
+!
+!-- Temporary solution to add LSM and radiation time series to the default
+!-- output
+    IF ( land_surface  .OR.  radiation )  THEN
+       IF ( TRIM( radiation_scheme ) == 'rrtmg' )  THEN
+          dots_num = dots_num + 15
+       ELSE
+          dots_num = dots_num + 11
+       ENDIF
+    ENDIF
+    
 !
 !-- If required, initialize urban surface model
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2269)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2270)
@@ -20,5 +20,4 @@
 ! Current revisions:
 ! -----------------
-! 
 ! 
 ! 
@@ -26,4 +25,13 @@
 ! -----------------
 ! $Id$
+! Revised parameterization of heat conductivity between skin layer and soil.
+! Temperature and moisture are now defined at the center of the layers.
+! Renamed veg_type to vegetation_type and pave_type to pavement_type_name
+! Renamed and reduced the number of look-up tables (vegetation_pars, soil_pars)
+! Revised land surface model initialization
+! Removed output of shf_eb and qsws_eb and removed _eb throughout code
+! Removed Clapp & Hornberger parameterization
+! 
+! 2249 2017-06-06 13:58:01Z sward
 ! 
 ! 2248 2017-06-06 13:52:54Z sward $
@@ -50,5 +58,5 @@
 ! 
 ! 2149 2017-02-09 16:57:03Z scharf
-! Land surface parameters II corrected for veg_type 18 and 19
+! Land surface parameters II corrected for vegetation_type 18 and 19
 ! 
 ! 2031 2016-10-21 15:11:58Z knoop
@@ -177,5 +185,5 @@
 ! ------------
 !> Land surface model, consisting of a solver for the energy balance at the
-!> surface and a four layer soil scheme. The scheme is similar to the TESSEL
+!> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
 !> scheme implemented in the ECMWF IFS model, with modifications according to
 !> H-TESSEL. The implementation is based on the formulation implemented in the 
@@ -195,5 +203,9 @@
 !>       with considerable precipitation.
 !> @todo Fix crashes with radiation_scheme == 'constant'
-!>
+!> @todo Revise calculation of f2 when wilting point is non-constant in the
+!>       soil
+!> @todo Solve problems with soil heat flux paramterization
+!> @todo Make pavement_depth variable for each surface element
+!> @todo Allow for zero soil moisture (currently, it is set to wilting point)
 !> @note No time step criterion is required as long as the soil layers do not
 !>       become too thin.
@@ -211,6 +223,6 @@
                initializing_actions, intermediate_timestep_count_max,          &
                land_surface, max_masks, precipitation, pt_surface,             &
-               rho_surface, roughness_length, surface_pressure,                &
-               timestep_scheme, tsc, z0h_factor, time_since_reference_point
+               rho_surface, surface_pressure, timestep_scheme, tsc,            &
+               time_since_reference_point
 
     USE indices,                                                               &
@@ -223,10 +235,11 @@
     USE radiation_model_mod,                                                   &
         ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
-               rad_lw_out_change_0, radiation_scheme, unscheduled_radiation_calls
+               rad_lw_out_change_0, radiation_scheme,                          &
+               unscheduled_radiation_calls
         
     USE statistics,                                                            &
         ONLY:  hom, statistic_regions
 
-    USE surface_mod,                                                        &
+    USE surface_mod,                                                           &
         ONLY :  surf_lsm_h, surf_lsm_v, surf_type
 
@@ -243,10 +256,10 @@
     REAL(wp), PARAMETER  ::                    &
               b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
-              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
-              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
-              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
+              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil (W/m/K)  
+              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix (W/m/K)
+              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water (W/m/K)
               psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
-              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
-              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
+              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil (J/m3/K)
+              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water (J/m3/K)
               m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
 
@@ -260,12 +273,22 @@
 !
 !-- LSM variables
-    INTEGER(iwp) :: nzb_soil = 0,  & !< bottom of the soil model (Earth's surface)
-                    nzt_soil = 0,  & !< top of the soil model
-                    nzs = 0,       & !< number of soil layers
-                    veg_type  = 2, & !< default NAMELIST veg_type_2d
-                    soil_type = 3    !< default NAMELIST soil_type_2d
+    CHARACTER(10) :: surface_type = 'netcdf'  !< general classification. Allowed are:
+                                             !< 'vegetation', 'pavement', ('building'), 
+                                             !< 'water', and 'netcdf'
+
+
+
+    INTEGER(iwp) :: nzb_soil = 0,        & !< bottom of the soil model (Earth's surface)
+                    nzt_soil = 7,        & !< top of the soil model
+                    nzs = 8,             & !< number of soil layers
+                    pavement_type = 1,   & !< default NAMELIST pavement_type                 
+                    soil_type = 3,       & !< default NAMELIST soil_type
+                    vegetation_type = 2, & !< default NAMELIST vegetation_type
+                    water_type = 1         !< default NAMELISt water_type
+                    
    
         
     LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
+               constant_roughness = .FALSE.,     & !< use fixed/dynamic roughness lengths for water surfaces
                force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
                aero_resist_kray = .TRUE.           !< flag to control parametrization of aerodynamic resistance at vertical surface elements
@@ -275,5 +298,5 @@
     REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
                 canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
-                c_surface   = 20000.0_wp,               & !< Surface (skin) heat capacity
+                c_surface = 9999999.9_wp,               & !< Surface (skin) heat capacity (J/m2/K)
                 drho_l_lv,                              & !< (rho_l * l_v)**-1
                 exn,                                    & !< value of the Exner function
@@ -283,7 +306,7 @@
                 hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
                 ke = 0.0_wp,                            & !< Kersten number
-                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil
-                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
-                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
+                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil (W/m/K)
+                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s (W/m2/K)
+                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u (W/m2/K)
                 leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
                 l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
@@ -292,7 +315,7 @@
                 m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
                 n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
-                pave_depth = 9999999.9_wp,              & !< depth of the pavement
-                pave_heat_capacity = 1.94E6_wp,         & !< volumetric heat capacity of pavement (e.g. roads)
-                pave_heat_conductivity = 1.00_wp,       & !< heat conductivity for pavements (e.g. roads)
+                pavement_depth = 9999999.9_wp,          & !< depth of the pavement
+                pavement_heat_capacity = 9999999.9_wp,  & !< volumetric heat capacity of pavement (e.g. roads) (J/m3/K)
+                pavement_heat_conduct  = 9999999.9_wp,  & !< heat conductivity for pavements (e.g. roads) (W/m/K)
                 q_s = 0.0_wp,                           & !< saturation specific humidity
                 residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
@@ -303,22 +326,31 @@
                 skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
                 vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
+                water_temperature = 9999999.9_wp,       & !< water temperature
                 wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
-                z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
-                z0h_eb = 9999999.9_wp,                  & !< NAMELIST z0h (lsm_par)
-                z0q_eb = 9999999.9_wp                     !< NAMELIST z0q (lsm_par)
-
+                z0_vegetation  = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
+                z0h_vegetation = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
+                z0q_vegetation = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
+                z0_pavement    = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
+                z0h_pavement   = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
+                z0q_pavement   = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
+                z0_water       = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
+                z0h_water      = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
+                z0q_water      = 9999999.9_wp             !< NAMELIST z0q (lsm_par)  
+                
+                
     REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil,         & !< 1/dz_soil
-                                            ddz_soil_stag,    & !< 1/dz_soil_stag
-                                            dz_soil,          & !< soil grid spacing (edge-edge)
-                                            dz_soil_stag,     & !< soil grid spacing (center-center)
+                                            ddz_soil_layer,   & !< 1/dz_soil_layer
+                                            dz_soil,          & !< soil grid spacing (center-center)
+                                            dz_soil_layer,    & !< soil grid spacing (edge-edge)
                                             root_extr           !< root extraction
 
 
                                             
-    REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp,               & !< distribution of root surface area to the individual soil layers
+    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
-                                   zs = 9999999.9_wp                             !< soil layer depths
-                                   
+                                   zs = 9999999.9_wp,                          & !< (NAMELIST) soil layer depths (center)
+                                   zs_layer = 9999999.9_wp                       !< soil layer depths (edge)
+                                 
 #if defined( __nopointer )
     TYPE(surf_type_lsm), TARGET  ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
@@ -359,40 +391,40 @@
     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_eb_h,     & !< liquid water reservoir (m), horizontal surface elements 
-                                      m_liq_eb_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   ::  &
                                       t_surface_v,    & !< surface temperature (K), vertical surface elements 
                                       t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements 
-                                      m_liq_eb_v,     & !< liquid water reservoir (m), vertical surface elements 
-                                      m_liq_eb_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_eb_h,     & !< liquid water reservoir (m), horizontal surface elements 
-                                      m_liq_eb_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_eb_h_1,   & !<
-                                      m_liq_eb_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_eb_v,     & !< liquid water reservoir (m), vertical surface elements 
-                                      m_liq_eb_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_eb_v_1,   & !<
-                                      m_liq_eb_v_2      !<
+                                      m_liq_v_1,   & !<
+                                      m_liq_v_2      !<
 #endif
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
 #else
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
 #endif
 
@@ -405,10 +437,10 @@
 #endif
 
-    TYPE(surf_type_lsm), TARGET ::  tm_liq_eb_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_eb_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 
@@ -421,44 +453,40 @@
               c_soil_av,        & !< average of c_soil
               c_veg_av,         & !< average of c_veg
-              ghf_eb_av,        & !< average of ghf_eb
+              ghf_av,        & !< average of ghf
               lai_av,           & !< average of lai
-              qsws_eb_av,       & !< average of qsws_eb
-              qsws_liq_eb_av,   & !< average of qsws_liq_eb
-              qsws_soil_eb_av,  & !< average of qsws_soil_eb
-              qsws_veg_eb_av,   & !< average of qsws_veg_eb
+              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
-              shf_eb_av           !< average of shf_eb
-
-
-!
-!-- Predefined Land surface classes (veg_type)
-    CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
-                                   'user defined              ',    & ! 0 
-                                   'crops, mixed farming      ',    & !  1
-                                   'short grass               ',    & !  2
-                                   'evergreen needleleaf trees',    & !  3
-                                   'deciduous needleleaf trees',    & !  4
-                                   'evergreen broadleaf trees ',    & !  5
-                                   'deciduous broadleaf trees ',    & !  6
-                                   'tall grass                ',    & !  7
-                                   'desert                    ',    & !  8
-                                   'tundra                    ',    & !  9
-                                   'irrigated crops           ',    & ! 10
-                                   'semidesert                ',    & ! 11
-                                   'ice caps and glaciers     ',    & ! 12
-                                   'bogs and marshes          ',    & ! 13
-                                   'inland water              ',    & ! 14
-                                   'ocean                     ',    & ! 15
-                                   'evergreen shrubs          ',    & ! 16
-                                   'deciduous shrubs          ',    & ! 17
-                                   'mixed forest/woodland     ',    & ! 18
-                                   'interrupted forest        ',    & ! 19
-                                   'pavements/roads           '     & ! 20
+              r_s_av              !< average of r_s
+                   
+
+!
+!-- Predefined Land surface classes (vegetation_type)
+    CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
+                                   'user defined              ',           & !  0 
+                                   'bare soil                 ',           & !  1                           
+                                   '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
                                                                  /)
 
 !
 !-- Soil model classes (soil_type)
-    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
+    CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
                                    'user defined',                  & ! 0 
                                    'coarse      ',                  & ! 1
@@ -467,131 +495,99 @@
                                    'fine        ',                  & ! 4
                                    'very fine   ',                  & ! 5
-                                   'organic     ',                  & ! 6
-                                   'loamy (CH)  '                   & ! 7
+                                   'organic     '                   & ! 6
                                                                  /)
-!
-!-- Land surface parameters according to the respective classes (veg_type)
-
-!
-!-- Land surface parameters I
-!--                          r_canopy_min,     lai,   c_veg,     g_d
-    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
-                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
-                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
-                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
-                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
-                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
-                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
-                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
-                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
-                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
-                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
-                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
-                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
-                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
-                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
-                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
-                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
-                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
-                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
-                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,  & ! 19
-                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp   & ! 20
-                                 /), (/ 4, 20 /) )
-
-!
-!-- Land surface parameters II          z0,         z0h,         z0q
-    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & 
-                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & !  1
-                                   0.20_wp,  0.20E-2_wp,  0.20E-2_wp,       & !  2
-                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  3
-                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  4
-                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  5
-                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  6
-                                   0.47_wp,  0.47E-2_wp,  0.47E-2_wp,       & !  7
-                                  0.013_wp, 0.013E-2_wp, 0.013E-2_wp,       & !  8
-                                  0.034_wp, 0.034E-2_wp, 0.034E-2_wp,       & !  9
-                                    0.5_wp,  0.50E-2_wp,  0.50E-2_wp,       & ! 10
-                                   0.17_wp,  0.17E-2_wp,  0.17E-2_wp,       & ! 11
-                                 1.3E-3_wp,   1.3E-4_wp,   1.3E-4_wp,       & ! 12
-                                   0.83_wp,  0.83E-2_wp,  0.83E-2_wp,       & ! 13
-                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 14
-                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 15
-                                   0.10_wp,  0.10E-2_wp,  0.10E-2_wp,       & ! 16
-                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & ! 17
-                                   2.00_wp,     2.00_wp,     2.00_wp,       & ! 18
-                                   1.10_wp,     1.10_wp,     1.10_wp,       & ! 19
-                                 1.0E-4_wp,   1.0E-5_wp,   1.0E-5_wp        & ! 20
-                                 /), (/ 3, 20 /) )
-
-!
-!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
-    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
-                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
-                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
-                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
-                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
-                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10 
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
-                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13 
-                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 14
-                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 15
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
-                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
-                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
-                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 19
-                                       0.0_wp,        0.0_wp, 0.00_wp      & ! 20
-                                      /), (/ 3, 20 /) )
-
-!
-!-- Root distribution (sum = 1)  level 1-8
-!--                               
-    REAL(wp), DIMENSION(0:7,1:20), PARAMETER :: root_distribution = RESHAPE( (/&
+
+!
+!-- 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
+                                                                 /)                                                                 
+                                                                 
+!
+!-- Water classes
+    CHARACTER(12), DIMENSION(0:4), PARAMETER :: water_type_name = (/ &
+                                   'user defined',                   & ! 0 
+                                   'lake        ',                   & ! 1
+                                   'river       ',                   & ! 2
+                                   'ocean       ',                   & ! 3
+                                   'pond        '                    & ! 4
+                                                                  /)                                                                  
+                                                                 
+                                                                 !
+!-- Land surface parameters according to the respective classes (vegetation_type)
+
+!
+!-- Land surface parameters
+!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface
+    REAL(wp), DIMENSION(0:9,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
+          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.5E-2_wp,   0.5E-4_wp,     0.0_wp,     0.0_wp, 0.00_wp, 0.00_wp, & !  1
+        180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  2
+        110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,   0.20_wp,  0.20E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  3
+        500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  4
+        500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  5
+        175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  6
+        240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  7
+        100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,   0.47_wp,  0.47E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  8
+        250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  0.013_wp, 0.013E-2_wp,   196.0_wp,   196.0_wp, 0.00_wp, 0.00_wp, & !  9
+         80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  0.034_wp, 0.034E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 10
+        180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,    0.5_wp,  0.50E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 11
+        150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,   0.17_wp,  0.17E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 12
+          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp,   1.3E-4_wp,   812.0_wp,   812.0_wp, 0.00_wp, 0.00_wp, & ! 13
+        240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,   0.83_wp,  0.83E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 14
+        225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,   0.10_wp,  0.10E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 15
+        225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 16
+        250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & ! 17
+        175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,   1.10_wp,     1.10_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp  & ! 18
+                                                               /), (/ 10, 18 /) )
+
+                                    
+!
+!-- Root distribution for default soil layer configuration (sum = 1)
+!--                                level 1 - level 8
+    REAL(wp), DIMENSION(0:7,1:18), PARAMETER :: root_distribution = RESHAPE( (/&
+                                   1.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
+                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 1
                                    0.035_wp, 0.069_wp, 0.069_wp, 0.108_wp,     &
-                                   0.195_wp, 0.214_wp, 0.284_wp, 0.026_wp,     & ! 1
+                                   0.195_wp, 0.214_wp, 0.284_wp, 0.026_wp,     & ! 2
                                    0.050_wp, 0.100_wp, 0.100_wp, 0.136_wp,     &
-                                   0.181_wp, 0.192_wp, 0.215_wp, 0.026_wp,     & ! 2
+                                   0.181_wp, 0.192_wp, 0.215_wp, 0.026_wp,     & ! 3
                                    0.038_wp, 0.075_wp, 0.075_wp, 0.111_wp,     &
-                                   0.185_wp, 0.203_wp, 0.273_wp, 0.040_wp,     & ! 3
+                                   0.185_wp, 0.203_wp, 0.273_wp, 0.040_wp,     & ! 4
                                    0.038_wp, 0.075_wp, 0.075_wp, 0.110_wp,     &
-                                   0.180_wp, 0.199_wp, 0.277_wp, 0.046_wp,     & ! 4
+                                   0.180_wp, 0.199_wp, 0.277_wp, 0.046_wp,     & ! 5
                                    0.035_wp, 0.069_wp, 0.069_wp, 0.105_wp,     &
-                                   0.180_wp, 0.201_wp, 0.295_wp, 0.046_wp,     & ! 5
+                                   0.180_wp, 0.201_wp, 0.295_wp, 0.046_wp,     & ! 6
                                    0.035_wp, 0.072_wp, 0.072_wp, 0.105_wp,     &
-                                   0.161_wp, 0.180_wp, 0.282_wp, 0.093_wp,     & ! 6
+                                   0.161_wp, 0.180_wp, 0.282_wp, 0.093_wp,     & ! 7
                                    0.040_wp, 0.077_wp, 0.077_wp, 0.112_wp,     &
-                                   0.176_wp, 0.192_wp, 0.266_wp, 0.060_wp,     & ! 7
+                                   0.176_wp, 0.192_wp, 0.266_wp, 0.060_wp,     & ! 8
                                    0.142_wp, 0.286_wp, 0.286_wp, 0.286_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 8
+                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 9
                                    0.068_wp, 0.134_wp, 0.134_wp, 0.177_wp,     &
-                                   0.214_wp, 0.203_wp, 0.070_wp, 0.000_wp,     & ! 9
+                                   0.214_wp, 0.203_wp, 0.070_wp, 0.000_wp,     & ! 10
                                    0.035_wp, 0.068_wp, 0.068_wp, 0.108_wp,     &
-                                   0.195_wp, 0.215_wp, 0.285_wp, 0.026_wp,     & ! 10
+                                   0.195_wp, 0.215_wp, 0.285_wp, 0.026_wp,     & ! 11
                                    0.025_wp, 0.048_wp, 0.048_wp, 0.078_wp,     &
-                                   0.147_wp, 0.175_wp, 0.353_wp, 0.126_wp,     & ! 11
+                                   0.147_wp, 0.175_wp, 0.353_wp, 0.126_wp,     & ! 12
                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 12
+                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 13
                                    0.036_wp, 0.072_wp, 0.072_wp, 0.103_wp,     &
-                                   0.163_wp, 0.180_wp, 0.273_wp, 0.074_wp,     & ! 13
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 14
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 15
+                                   0.163_wp, 0.180_wp, 0.273_wp, 0.074_wp,     & ! 14
+                                   0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
+                                   0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 15
                                    0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
                                    0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 16
-                                   0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
-                                   0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 17
                                    0.028_wp, 0.055_wp, 0.055_wp, 0.087_wp,     &
-                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp,     & ! 18
+                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp,     & ! 17
                                    0.028_wp, 0.055_wp, 0.055_wp, 0.087_wp,     &
-                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp,     & ! 19
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp      & ! 20
-                                 /), (/ 8, 20 /) )
+                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp      & ! 18
+                                 /), (/ 8, 18 /) )
 
 
@@ -600,28 +596,37 @@
 
 !
-!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
-    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
-                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
-                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
-                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
-                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
-                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
-                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
-                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
-                                 /), (/ 4, 7 /) )
-
-!
-!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res  
-    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
-                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
-                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
-                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
-                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
-                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
-                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
-                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
-                                 /), (/ 4, 7 /) )
-
-
+!-- 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 /) )
+
+!
+!-- TO BE FILLED
+!-- Pavement parameters  depth,        z0,       z0h, lambda_h,      rho_c  
+    REAL(wp), DIMENSION(0:4,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ &
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 1
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 2
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 3                                  
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 4
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 5
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 6
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp  & ! 7
+                                                                 /), (/ 5, 7 /) )                             
+
+!
+!-- TO BE FILLED
+!-- Water parameters                    temperature,     z0,      z0h
+    REAL(wp), DIMENSION(0:2,1:4), PARAMETER :: water_pars = RESHAPE( (/ &
+                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 1
+                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 2
+                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 3
+                                           283.0_wp, 10.01_wp, 0.001_wp & ! 4
+                                                                /), (/ 3, 4 /) )                                                                     
+                                                                                                                                      
     SAVE
 
@@ -637,5 +642,5 @@
            lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
            lsm_swap_timelevel, lsm_read_restart_data, lsm_last_actions
-!
+! !vegetat
 !-- Public parameters, constants and initial values
     PUBLIC aero_resist_kray, skip_time_do_lsm
@@ -754,7 +759,7 @@
           unit = 'K'   
              
-       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
-              'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
-              'r_a*', 'r_s*', 'shf_eb*' )
+       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: "' //         &
@@ -783,20 +788,15 @@
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
+          IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
              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_eb*'  .AND.  .NOT.  land_surface )  THEN
+          IF ( TRIM( var ) == 'm_liq*'  .AND.  .NOT.  land_surface )  THEN
              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_eb*'  .AND.  .NOT.  land_surface )  THEN
-             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_eb*'  .AND.  .NOT. land_surface )   &
+          IF ( TRIM( var ) == 'qsws_liq*'  .AND.  .NOT. land_surface )   &
           THEN
              message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -804,5 +804,5 @@
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
+          IF ( TRIM( var ) == 'qsws_soil*'  .AND.  .NOT.  land_surface ) &
           THEN
              message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -810,5 +810,5 @@
              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
+          IF ( TRIM( var ) == 'qsws_veg*'  .AND.  .NOT. land_surface )   &
           THEN
              message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -833,13 +833,11 @@
           IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
           IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
-          IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
-          IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
-          IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
-          IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
-          IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
-          IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
+          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 ) == 'shf_eb*')  unit = 'W/m2'
              
        CASE DEFAULT
@@ -888,8 +886,8 @@
              dopr_index(var_count) = 89
              dopr_unit     = 'K'
-             hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+             hom(0:nzs-1,2,89,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
              IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
                 dopr_initial_index(var_count) = 90
-                hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                hom(0:nzs-1,2,90,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
                 data_output_pr(var_count)     = data_output_pr(var_count)(2:)
              ENDIF
@@ -906,8 +904,8 @@
              dopr_index(var_count) = 91
              dopr_unit     = 'm3/m3'
-             hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+             hom(0:nzs-1,2,91,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
              IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
                 dopr_initial_index(var_count) = 92
-                hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                hom(0:nzs-1,2,92,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
                 data_output_pr(var_count)     = data_output_pr(var_count)(2:)
              ENDIF
@@ -943,4 +941,17 @@
 
     INTEGER(iwp) ::  k        !< running index, z-dimension
+
+!
+!-- Check for a valid setting of surface_type. The default value is 'netcdf'.
+!-- In that case, the surface types are read from NetCDF file
+    IF ( TRIM( surface_type ) /= 'vegetation'  .AND.                           &
+         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 ) // '"'
+                CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
+    ENDIF
+
 !
 !-- Dirichlet boundary conditions are required as the surface fluxes are
@@ -948,6 +959,6 @@
 !-- model
     IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
-       message_string = 'lsm requires setting of'//                         &
-                        'bc_pt_b = "dirichlet" and '//                      &
+       message_string = 'lsm requires setting of'//                            &
+                        'bc_pt_b = "dirichlet" and '//                         &
                         'bc_q_b  = "dirichlet"'
        CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
@@ -955,84 +966,159 @@
 
     IF (  .NOT.  constant_flux_layer )  THEN
-       message_string = 'lsm requires '//                                   &
+       message_string = 'lsm requires '//                                      &
                         'constant_flux_layer = .T.'
        CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
     ENDIF
 
-    IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
-           most_method == 'lookup' )  THEN
-        WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
-                                   'allowed in combination with ',          &
-                                   'most_method = ', most_method
-       CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
+    IF ( TRIM( surface_type ) == 'vegetation' )  THEN
+    
+       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'//       &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( leaf_area_index == 9999999.9_wp )  THEN
+             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 )
+          ENDIF
+
+          IF ( vegetation_coverage == 9999999.9_wp )  THEN
+             message_string = 'vegetation_type = 0 (user_defined)'//              &
+                              'requires setting of vegetation_coverage'//         &
+                              '/= 9999999.9'
+                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
+             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 )
+          ENDIF
+
+          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
+             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 )
+          ENDIF
+
+          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
+             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 )
+          ENDIF
+
+          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
+             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 )
+          ENDIF
+
+          IF ( z0_vegetation == 9999999.9_wp )  THEN
+             message_string = 'vegetation_type = 0 (user_defined)'//              &
+                              'requires setting of z0_vegetation'//               &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( z0h_vegetation == 9999999.9_wp )  THEN
+             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
+       
     ENDIF
-
-    IF ( veg_type == 0 )  THEN
-       IF ( min_canopy_resistance == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user defined)'//                  &
-                           'requires setting of min_canopy_resistance'//    &
-                           '/= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( leaf_area_index == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of leaf_area_index'//          &
-                           '/= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( vegetation_coverage == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of vegetation_coverage'//      &
-                           '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of'//                          &
-                           'canopy_resistance_coefficient /= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of lambda_surface_stable'//    &
-                           '/= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of lambda_surface_unstable'//  &
-                           '/= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of f_shortwave_incoming'//     &
-                           '/= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( z0_eb == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of z0_eb'//                    &
-                           '/= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( z0h_eb == 9999999.9_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of z0h_eb'//                   &
-                           '/= 9999999.9'
-          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-       ENDIF
-
-
+    
+    IF ( TRIM( surface_type ) == 'water' )  THEN
+
+       IF ( TRIM( most_method ) == 'lookup' )  THEN    
+          WRITE( message_string, * ) 'surface_type = ', surface_type,          &
+                                     ' is not allowed in combination with ',   &
+                                     'most_method = ', most_method
+          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( water_type == 0 )  THEN  
+       
+          IF ( z0_water == 9999999.9_wp )  THEN
+             message_string = 'water_type = 0 (user_defined)'//                &
+                              'requires setting of z0_water'//                 &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( z0h_water == 9999999.9_wp )  THEN
+             message_string = 'water_type = 0 (user_defined)'//                &
+                              'requires setting of z0h_water'//                &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+          IF ( water_temperature == 9999999.9_wp )  THEN
+             message_string = 'water_type = 0 (user_defined)'//                &
+                              'requires setting of water_temperature'//        &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
+          ENDIF       
+          
+       ENDIF
+       
     ENDIF
-
+    
+    
+     
+    IF ( TRIM( surface_type ) == 'pavement' )  THEN
+
+       IF ( pavement_type == 0 )  THEN  
+       
+          IF ( z0_pavement == 9999999.9_wp )  THEN
+             message_string = 'pavement_type = 0 (user_defined)'//             &
+                              'requires setting of z0_pavement'//              &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( z0h_pavement == 9999999.9_wp )  THEN
+             message_string = 'pavement_type = 0 (user_defined)'//             &
+                              'requires setting of z0h_pavement'//             &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
+          ENDIF
+
+           IF ( pavement_depth == 9999999.9_wp )  THEN
+             message_string = 'pavement_type = 0 (user_defined)'//             &
+                              'requires setting of pavement_depth'//           &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
+          ENDIF         
+          
+           IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
+             message_string = 'pavement_type = 0 (user_defined)'//             &
+                              'requires setting of pavement_heat_conduct'//    &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
+          ENDIF 
+          
+           IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
+             message_string = 'pavement_type = 0 (user_defined)'//             &
+                              'requires setting of pavement_heat_capacity'//   &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
+          ENDIF 
+
+       ENDIF 
+    
+    ENDIF
+    
     IF ( soil_type == 0 )  THEN
 
@@ -1094,22 +1180,4 @@
 
     ENDIF
-!
-!-- Check for proper setting of soil moisture, must not be larger than its 
-!-- saturation value. 
-    DO  k = nzb_soil, nzt_soil
-       IF ( soil_moisture(k) > m_soil_pars(0,soil_type) )  THEN 
-          message_string = 'soil_moisture must not exceed its saturation' // &
-                           ' value'
-          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-       ENDIF
-    ENDDO
-
-    IF (  .NOT.  radiation )  THEN
-       message_string = 'lsm requires '//                                   &
-                        'radiation = .T.'
-       CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
-    ENDIF
-       
-
 
 !
@@ -1117,5 +1185,5 @@
 !-- root fraction is prescribed
     nzb_soil = 0
-    nzt_soil = 0
+    nzt_soil = -1
     IF ( ALL( zs == 9999999.9_wp ) )  THEN
        nzt_soil = 7
@@ -1136,13 +1204,24 @@
     ENDIF
 
-    IF ( veg_type == 0 )  THEN
+    IF ( vegetation_type == 0 )  THEN
        IF ( SUM( root_fraction(nzb_soil:nzt_soil) ) /= 1.0_wp )  THEN
-          message_string = 'veg_type = 0 (user_defined)'//                  &
-                           'requires setting of root_fraction'//            &
+          message_string = 'vegetation_type = 0 (user_defined)'//              &
+                           'requires setting of root_fraction'//               &
                            '/= 9999999.9 and SUM(root_fraction) = 1'
           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
-    ENDIF
+    ENDIF    
     
+    
+!
+!-- Check for proper setting of soil moisture, must not be larger than its 
+!-- saturation value. 
+    DO  k = nzb_soil, nzt_soil
+       IF ( soil_moisture(k) > soil_pars(4,soil_type) )  THEN 
+          message_string = 'soil_moisture must not exceed its saturation' // &
+                           ' value'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+    ENDDO
   
  END SUBROUTINE lsm_check_parameters
@@ -1155,4 +1234,6 @@
  SUBROUTINE lsm_energy_balance( horizontal, l )
 
+ 
+    USE pegrid
 
     IMPLICIT NONE
@@ -1183,11 +1264,11 @@
                 coef_1,      & !< coef. for prognostic equation
                 coef_2,      & !< coef. for prognostic equation
-                f_qsws,      & !< factor for qsws_eb
-                f_qsws_veg,  & !< factor for qsws_veg_eb
-                f_qsws_soil, & !< factor for qsws_soil_eb
-                f_qsws_liq,  & !< factor for qsws_liq_eb
-                f_shf,       & !< factor for shf_eb
-                lambda_surface, & !< Current value of lambda_surface
-                m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
+                f_qsws,      & !< factor for qsws
+                f_qsws_veg,  & !< factor for qsws_veg
+                f_qsws_soil, & !< factor for qsws_soil
+                f_qsws_liq,  & !< factor for qsws_liq
+                f_shf,       & !< factor for shf
+                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
@@ -1196,7 +1277,7 @@
     TYPE(surf_type_lsm), POINTER ::  surf_t_surface_p
     TYPE(surf_type_lsm), POINTER ::  surf_tt_surface_m
-    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_eb
-    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_eb_p
-    TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_eb_m
+    TYPE(surf_type_lsm), POINTER ::  surf_m_liq
+    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_p
+    TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_m
 
     TYPE(surf_type_lsm), POINTER ::  surf_m_soil
@@ -1211,7 +1292,7 @@
        surf_t_surface_p  => t_surface_h_p
        surf_tt_surface_m => tt_surface_h_m
-       surf_m_liq_eb     => m_liq_eb_h
-       surf_m_liq_eb_p   => m_liq_eb_h_p
-       surf_tm_liq_eb_m  => tm_liq_eb_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
@@ -1226,7 +1307,7 @@
        surf_t_surface_p  => t_surface_v_p(l)
        surf_tt_surface_m => tt_surface_v_m(l)
-       surf_m_liq_eb     => m_liq_eb_v(l)
-       surf_m_liq_eb_p   => m_liq_eb_v_p(l)
-       surf_tm_liq_eb_m  => tm_liq_eb_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)
@@ -1263,20 +1344,41 @@
 !--    radiation can be obtained at respective height level
        k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )
-!
-!--    Set lambda_surface according to stratification between skin layer and soil
-       IF (  .NOT.  surf%pave_surface(m) )  THEN
-
-          c_surface_tmp = c_surface
-
-          IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))  THEN
-             lambda_surface = surf%lambda_surface_s(m)
-          ELSE
-             lambda_surface = surf%lambda_surface_u(m)
-          ENDIF
+
+!
+!--    Define heat conductivity between surface and soil depending on surface
+!--    type. For vegetation, a skin layer parameterization is used. For bare
+!--    soil and pavements, no skin layer is applied. In these cases, the
+!--    temperature is assumed to be constant between the surface and the first
+!--    soil layer. The heat conductivity is then derived from the soil/
+!--    pavement properties. For water surfaces, the conductivity is already set
+!--    to 1E10.
+!--    Moreover, the heat capacity is set. For bare soil the heat capacity is
+!--    the capacity of the uppermost soil layer
+       IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN 
+       
+          lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
+                         lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
+                         
+          ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) /   &
+                                                     surf%m_sat(nzb_soil,m) ) )                   
+                         
+          lambda_surface = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )&
+                           * ddz_soil_layer(nzb_soil) * 2.0_wp
+          
+          c_surface_tmp = (rho_c_soil * ( 1.0_wp - surf%m_sat(nzb_soil,m) )    &
+                           + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) )    &
+                           * dz_soil_layer(nzb_soil) * 0.25_wp               
+
+       ELSEIF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m)) &
+       THEN
+
+          lambda_surface = surf%lambda_surface_s(m)
+          c_surface_tmp = surf%c_surface(m) 
+          
        ELSE
 
-          c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
-          lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
-
+          lambda_surface = surf%lambda_surface_u(m)
+          c_surface_tmp = surf%c_surface(m) 
+          
        ENDIF
 
@@ -1333,12 +1435,16 @@
        DO  ks = nzb_soil, nzt_soil
            m_total = m_total + surf%root_fr(ks,m)                              &
-                     * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(m) )
+                     * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
        ENDDO 
 
-       IF ( m_total > surf%m_wilt(m)  .AND.                                    &
-            m_total < surf%m_fc(m) )  THEN
-          f2 = ( m_total - surf%m_wilt(m) ) /                                  &
-               ( surf%m_fc(m) - surf%m_wilt(m) )
-       ELSEIF ( m_total >= surf%m_fc(m) )  THEN
+!
+!--    The calculation of f2 is based on only one wilting point value for all
+!--    soil layers. The value at k=nzb_soil is used here as a proxy but might
+!--    need refinement in the future.
+       IF ( m_total > surf%m_wilt(nzb_soil,m)  .AND.                           &
+            m_total < surf%m_fc(nzb_soil,m) )  THEN
+          f2 = ( m_total - surf%m_wilt(nzb_soil,m) ) /                         &
+               ( surf%m_fc(nzb_soil,m) - surf%m_wilt(nzb_soil,m) )
+       ELSEIF ( m_total >= surf%m_fc(nzb_soil,m) )  THEN
           f2 = 1.0_wp
        ELSE
@@ -1366,19 +1472,15 @@
                               ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
 !
-!--    Third step: calculate bare soil resistance r_soil. The Clapp & 
-!--    Hornberger parametrization does not consider c_veg.
-       IF ( surf%soil_type_2d(m) /= 7 )  THEN
-          m_min = surf%c_veg(m) * surf%m_wilt(m) +                             &
-                         ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(m)
-       ELSE
-          m_min = surf%m_wilt(m)
-       ENDIF
-
-       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) / ( surf%m_fc(m) - m_min )
+!--    Third step: calculate bare soil resistance r_soil.
+       m_min = surf%c_veg(m) * surf%m_wilt(nzb_soil,m) +                       &
+                         ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(nzb_soil,m)
+
+
+       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     )
 
        surf%r_soil(m) = surf%r_soil_min(m) / f2
-
+       
 !
 !--    Calculate the maximum possible liquid water amount on plants and
@@ -1388,11 +1490,11 @@
 !--    Noilhan & Planton (1989), while the ECMWF formulation is used for
 !--    vegetated surfaces and bare soils.
-       IF ( surf%pave_surface(m) )  THEN
-          m_liq_eb_max = m_max_depth * 5.0_wp
-          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq_eb%var_1d(m) / m_liq_eb_max)**0.67 )
+       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 )
        ELSE
-          m_liq_eb_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)&
-                            + ( 1.0_wp - surf%c_veg(m) ) )
-          surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq_eb%var_1d(m) / m_liq_eb_max )
+          m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)              &
+                      + ( 1.0_wp - surf%c_veg(m) ) )
+          surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq%var_1d(m) / m_liq_max )
        ENDIF
 !
@@ -1415,5 +1517,5 @@
           f_qsws_soil = 0.0_wp
           f_qsws_liq  = rho_lv / surf%r_a(m)
-       ELSEIF ( surf%pave_surface (m) )  THEN
+       ELSEIF ( surf%pavement_surface (m) )  THEN
           f_qsws_veg  = 0.0_wp
           f_qsws_soil = 0.0_wp
@@ -1428,10 +1530,4 @@
                               surf%r_a(m)
        ENDIF
-!
-!--       If soil moisture is below wilting point, plants do no longer
-!--       transpirate.
-!           IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
-!              f_qsws_veg = 0.0_wp
-!           ENDIF
 
        f_shf  = rho_cp / surf%r_a(m)
@@ -1535,27 +1631,24 @@
                      ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
 
-       surf%ghf_eb(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) )
 
-       surf%shf_eb(m) = - f_shf * ( pt1 - surf%pt_surface(m) )
-
-       surf%shf(m) = surf%shf_eb(m) / rho_cp
+       surf%shf(m) = - f_shf * ( pt1 - surf%pt_surface(m) ) / cp
+
 
        IF ( humidity )  THEN
-          surf%qsws_eb(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(m) = surf%qsws_eb(m) / rho_lv
-
-          surf%qsws_veg_eb(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_eb(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_eb(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) )
@@ -1564,5 +1657,5 @@
 !
 !--    Calculate the true surface resistance
-       IF ( surf%qsws_eb(m) == 0.0_wp )  THEN
+       IF ( .NOT.  humidity )  THEN
           surf%r_s(m) = 1.0E10_wp
        ELSE
@@ -1570,5 +1663,5 @@
                           * surf_t_surface%var_1d(m) - dq_s_dt *               &
                             surf_t_surface_p%var_1d(m) ) /                     &
-                            surf%qsws_eb(m) - surf%r_a(m)
+                            surf%qsws(m) - surf%r_a(m)
        ENDIF
 
@@ -1578,6 +1671,6 @@
        IF ( humidity )  THEN
 !
-!--       If precipitation is activated, add rain water to qsws_liq_eb
-!--       and qsws_soil_eb according the the vegetation coverage.
+!--       If precipitation is activated, add rain water to qsws_liq
+!--       and qsws_soil according the the vegetation coverage.
 !--       precipitation_rate is given in mm.
           IF ( precipitation )  THEN
@@ -1587,12 +1680,12 @@
 !--          Otherwise, add the water to soil. In case of
 !--          pavements, the exceeding water amount is implicitely removed 
-!--          as runoff as qsws_soil_eb is then not used in the soil model
-             IF ( surf_m_liq_eb%var_1d(m) /= m_liq_eb_max )  THEN
-                surf%qsws_liq_eb(m) = surf%qsws_liq_eb(m)                      &
+!--          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%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_eb(m) = surf%qsws_soil_eb(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)                              &
@@ -1602,5 +1695,5 @@
 !--          Add precipitation to bare soil according to the bare soil
 !--          coverage.
-             surf%qsws_soil_eb(m) = surf%qsws_soil_eb(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)                                &
@@ -1610,48 +1703,50 @@
 !
 !--       If the air is saturated, check the reservoir water level
-          IF ( surf%qsws_eb(m) < 0.0_wp )  THEN
-!
-!--          Check if reservoir is full (avoid values > m_liq_eb_max)
-!--          In that case, qsws_liq_eb goes to qsws_soil_eb. In this 
-!--          case qsws_veg_eb is zero anyway (because c_liq = 1),       
+          IF ( surf%qsws(m) < 0.0_wp )  THEN
+!
+!--          Check if reservoir is full (avoid values > m_liq_max)
+!--          In that case, qsws_liq goes to qsws_soil. In this 
+!--          case qsws_veg is zero anyway (because c_liq = 1),       
 !--          so that tend is zero and no further check is needed
-             IF ( surf_m_liq_eb%var_1d(m) == m_liq_eb_max )  THEN
-                surf%qsws_soil_eb(m) = surf%qsws_soil_eb(m) + surf%qsws_liq_eb(m)
-
-                surf%qsws_liq_eb(m)  = 0.0_wp
+             IF ( surf_m_liq%var_1d(m) == m_liq_max )  THEN
+                surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
+
+                surf%qsws_liq(m)  = 0.0_wp
              ENDIF
 
 !
-!--          In case qsws_veg_eb becomes negative (unphysical behavior), 
+!--          In case qsws_veg becomes negative (unphysical behavior), 
 !--          let the water enter the liquid water reservoir as dew on the
 !--          plant
-             IF ( surf%qsws_veg_eb(m) < 0.0_wp )  THEN
-                surf%qsws_liq_eb(m) = surf%qsws_liq_eb(m) + surf%qsws_veg_eb(m)
-                surf%qsws_veg_eb(m) = 0.0_wp
+             IF ( surf%qsws_veg(m) < 0.0_wp )  THEN
+                surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
+                surf%qsws_veg(m) = 0.0_wp
              ENDIF
           ENDIF                    
  
-          tend = - surf%qsws_liq_eb(m) * drho_l_lv
-          surf_m_liq_eb_p%var_1d(m) = surf_m_liq_eb%var_1d(m) + dt_3d *        &
+          surf%qsws(m) = surf%qsws(m) / l_v
+ 
+          tend = - surf%qsws_liq(m) * drho_l_lv
+          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *        &
                                         ( tsc(2) * tend +                      &
-                                          tsc(3) * surf_tm_liq_eb_m%var_1d(m) )
+                                          tsc(3) * surf_tm_liq_m%var_1d(m) )
 !
 !--       Check if reservoir is overfull -> reduce to maximum
 !--       (conservation of water is violated here)
-          surf_m_liq_eb_p%var_1d(m) = MIN( surf_m_liq_eb_p%var_1d(m),m_liq_eb_max )
+          surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )
 
 !
 !--       Check if reservoir is empty (avoid values < 0.0)
 !--       (conservation of water is violated here)
-          surf_m_liq_eb_p%var_1d(m) = MAX( surf_m_liq_eb_p%var_1d(m), 0.0_wp )
-!
-!--       Calculate m_liq_eb tendencies for the next Runge-Kutta step
+          surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
+!
+!--       Calculate m_liq tendencies for the next Runge-Kutta step
           IF ( timestep_scheme(1:5) == 'runge' )  THEN
              IF ( intermediate_timestep_count == 1 )  THEN
-                surf_tm_liq_eb_m%var_1d(m) = tend
+                surf_tm_liq_m%var_1d(m) = tend
              ELSEIF ( intermediate_timestep_count <                            &
                       intermediate_timestep_count_max )  THEN
-                surf_tm_liq_eb_m%var_1d(m) = -9.5625_wp * tend +               &
-                                              5.3125_wp * surf_tm_liq_eb_m%var_1d(m)
+                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +               &
+                                              5.3125_wp * surf_tm_liq_m%var_1d(m)
              ENDIF
           ENDIF
@@ -1684,5 +1779,5 @@
 !
 !-- Calculate new roughness lengths (for water surfaces only)
-    IF ( horizontal )  CALL calc_z0_water_surface
+    IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
 
     CONTAINS
@@ -1799,5 +1894,5 @@
        ENDIF
 
-       WRITE( io, 4 ) TRIM( veg_type_name(veg_type) ),                         &
+       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 ),              &
@@ -1829,5 +1924,7 @@
     SUBROUTINE lsm_init
     
-
+       USE control_parameters,                                                 &
+           ONLY:  message_string
+    
        IMPLICIT NONE
 
@@ -1864,19 +1961,12 @@
        tt_soil_h_m%var_2d    = 0.0_wp
        tm_soil_h_m%var_2d    = 0.0_wp
-       tm_liq_eb_h_m%var_1d  = 0.0_wp
+       tm_liq_h_m%var_1d  = 0.0_wp
        surf_lsm_h%c_liq = 0.0_wp
 
-       surf_lsm_h%ghf_eb = 0.0_wp
-       surf_lsm_h%shf_eb = rho_cp * surf_lsm_h%shf
-
-       IF ( humidity )  THEN
-          surf_lsm_h%qsws_eb = rho_lv * surf_lsm_h%qsws
-       ELSE
-          surf_lsm_h%qsws_eb = 0.0_wp
-       ENDIF
-
-       surf_lsm_h%qsws_liq_eb  = 0.0_wp
-       surf_lsm_h%qsws_soil_eb = 0.0_wp
-       surf_lsm_h%qsws_veg_eb  = 0.0_wp
+       surf_lsm_h%ghf = 0.0_wp
+
+       surf_lsm_h%qsws_liq  = 0.0_wp
+       surf_lsm_h%qsws_soil = 0.0_wp
+       surf_lsm_h%qsws_veg  = 0.0_wp
 
        surf_lsm_h%r_a        = 50.0_wp
@@ -1890,19 +1980,12 @@
           tt_soil_v_m(l)%var_2d    = 0.0_wp
           tm_soil_v_m(l)%var_2d    = 0.0_wp
-          tm_liq_eb_v_m(l)%var_1d  = 0.0_wp
+          tm_liq_v_m(l)%var_1d  = 0.0_wp
           surf_lsm_v(l)%c_liq = 0.0_wp
 
-          surf_lsm_v(l)%ghf_eb = 0.0_wp
-          surf_lsm_v(l)%shf_eb = rho_cp * surf_lsm_v(l)%shf
-
-          IF ( humidity )  THEN
-             surf_lsm_v(l)%qsws_eb = rho_lv * surf_lsm_v(l)%qsws
-          ELSE
-             surf_lsm_v(l)%qsws_eb = 0.0_wp
-          ENDIF
-
-          surf_lsm_v(l)%qsws_liq_eb  = 0.0_wp
-          surf_lsm_v(l)%qsws_soil_eb = 0.0_wp
-          surf_lsm_v(l)%qsws_veg_eb  = 0.0_wp
+          surf_lsm_v(l)%ghf = 0.0_wp
+
+          surf_lsm_v(l)%qsws_liq  = 0.0_wp
+          surf_lsm_v(l)%qsws_soil = 0.0_wp
+          surf_lsm_v(l)%qsws_veg  = 0.0_wp
 
           surf_lsm_v(l)%r_a        = 50.0_wp
@@ -1911,14 +1994,21 @@
           surf_lsm_v(l)%r_soil     = 0.0_wp
        ENDDO
-
-        
+   
 !
 !--    Allocate 3D soil model arrays
 !--    First, for horizontal surfaces
+       ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
+       ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
+       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
+       ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
+       ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
+       ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
+       ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
+       ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns)      )
+       ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
+       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
        ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
-       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
-       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
-
-       surf_lsm_h%lambda_h = 0.0_wp
+    
+       surf_lsm_h%lambda_h     = 0.0_wp
 !
 !--    If required, allocate humidity-related variables for the soil model
@@ -1932,9 +2022,18 @@
 !--    For vertical surfaces
        DO  l = 0, 3
+          ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
+          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
+          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
+          ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
+          ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
+          ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
+          ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
+          ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)      )
+          ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
+          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )  
           ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
-          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
-          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
-
-          surf_lsm_v(l)%lambda_h = 0.0_wp
+
+          surf_lsm_v(l)%lambda_h     = 0.0_wp 
+          
 !
 !--       If required, allocate humidity-related variables for the soil model
@@ -1949,21 +2048,76 @@
 !
 !--    Calculate grid spacings. Temperature and moisture are defined at
-!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
-!--    at the centers
-       dz_soil(nzb_soil) = zs(nzb_soil)
-
+!--    the center of the soil layers, whereas gradients/fluxes are 
+!--    defined at the edges (_layer)
+       DO  k = nzb_soil, nzt_soil-1
+          dz_soil(k) = zs(k+1) - zs(k)
+          
+          IF ( dz_soil(k) == 0.0_wp )  THEN
+             message_string = 'invalid soil layer configuration found ' //     &
+                              '(dz_soil(k) = 0.0)'
+             CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+       ENDDO
+
+       dz_soil_layer(nzb_soil) = 2.0_wp * zs(nzb_soil)
+       zs_layer(nzb_soil)      = 2.0_wp * zs(nzb_soil)
+
+
+       
+       
        DO  k = nzb_soil+1, nzt_soil
-          dz_soil(k) = zs(k) - zs(k-1)
+          dz_soil_layer(k) = ( zs(k) - zs_layer(k-1) ) * 2.0_wp
+          zs_layer(k) = zs_layer(k-1) + dz_soil_layer(k)
        ENDDO
-       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
-
-       DO  k = nzb_soil, nzt_soil-1
-          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
+       dz_soil_layer(nzt_soil+1) = 0.5_wp * dz_soil_layer(nzt_soil)
+       dz_soil(nzt_soil) = zs_layer(k-1) + dz_soil_layer(k) - zs(nzt_soil)
+
+       
+       ddz_soil       = 1.0_wp / dz_soil
+       ddz_soil_layer = 1.0_wp / dz_soil_layer
+
+
+!
+!--    Set flag parameter for the prescribed surface type depending on user 
+!--    input.
+       DO  m = 1, surf_lsm_h%ns
+
+          SELECT CASE ( TRIM( surface_type ) )
+          
+             CASE ( 'vegetation' )
+         
+                surf_lsm_h%vegetation_surface(m)    = .TRUE.
+   
+             CASE ( 'water' )
+             
+                surf_lsm_h%water_surface(m)      = .TRUE.
+
+             CASE ( 'pavement' )
+             
+                surf_lsm_h%pavement_surface(m)   = .TRUE.
+    
+          END SELECT
+  
        ENDDO
-       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
-
-       ddz_soil      = 1.0_wp / dz_soil
-       ddz_soil_stag = 1.0_wp / dz_soil_stag
-
+       
+       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
@@ -1989,23 +2143,119 @@
 
           IF ( saturation_moisture == 9999999.9_wp )  THEN
-             saturation_moisture = m_soil_pars(0,soil_type)           
+             saturation_moisture = soil_pars(4,soil_type)           
           ENDIF
 
           IF ( field_capacity == 9999999.9_wp )  THEN
-             field_capacity = m_soil_pars(1,soil_type)           
+             field_capacity = soil_pars(5,soil_type)           
           ENDIF
 
           IF ( wilting_point == 9999999.9_wp )  THEN
-             wilting_point = m_soil_pars(2,soil_type)            
+             wilting_point = soil_pars(6,soil_type)            
           ENDIF
 
           IF ( residual_moisture == 9999999.9_wp )  THEN
-             residual_moisture = m_soil_pars(3,soil_type)       
+             residual_moisture = soil_pars(7,soil_type)       
           ENDIF
 
        ENDIF    
-
-!
-!--    Map values to the respective 2D arrays
+  
+!
+!--    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   
+          
+!
+!--       The heat conductivity between canopy and first soil layer must
+!--       depend on the depth of the soil layer. The value in the lookup
+!--       table refers to the first soil temperature at zs = -0.005 m. Lambda 
+!--       is thus interpolated to the actual soil layer depth used.
+          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
+             lambda_surface_stable = vegetation_pars(6,vegetation_type)        &
+                                     * 0.005_wp * ddz_soil_layer(nzb_soil)     &
+                                     * 2.0_wp    
+          ENDIF
+
+          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
+             lambda_surface_unstable = vegetation_pars(7,vegetation_type)      &
+                                       * 0.005_wp  * ddz_soil_layer(nzb_soil)  &
+                                       * 2.0_wp             
+          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
+              
+       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  
+
+       ENDIF 
+       
+       IF ( pavement_type /= 0 )  THEN  
+
+          IF ( pavement_depth == 9999999.9_wp )  THEN
+             pavement_depth = pavement_pars(0,pavement_type)
+          ELSE
+             pavement_depth = MAX( zs(nzb_soil), pavement_depth )
+          ENDIF
+
+          IF ( z0_pavement == 9999999.9_wp )  THEN
+             z0_pavement  = pavement_pars(1,pavement_type) 
+          ENDIF
+
+          IF ( z0h_pavement == 9999999.9_wp )  THEN
+             z0h_pavement = pavement_pars(2,pavement_type)
+          ENDIF
+
+          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
+             pavement_heat_conduct = pavement_pars(3,pavement_type)
+          ENDIF
+
+          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
+             pavement_heat_capacity = pavement_pars(4,pavement_type)
+          ENDIF   
+    
+       ENDIF
+!
+!--    Map values to the respective 2D/3D arrays
 !--    Horizontal surfaces
        surf_lsm_h%alpha_vg      = alpha_vangenuchten
@@ -2032,6 +2282,4 @@
        ENDDO
 
-
-
 !
 !--    Initial run actions
@@ -2039,32 +2287,79 @@
 !
 !--       First, for horizontal surfaces
-          t_soil_h%var_2d    = 0.0_wp
-          m_soil_h%var_2d    = 0.0_wp
-          m_liq_eb_h%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)
-          DO  k = nzb_soil, nzt_soil
-             t_soil_h%var_2d(k,:)      = soil_temperature(k)
-             m_soil_h%var_2d(k,:)      = MAX(soil_moisture(k),surf_lsm_h%m_wilt(:))
-             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
-          ENDDO
-          t_soil_h%var_2d(nzt_soil+1,:) = soil_temperature(nzt_soil+1)
-
-!
-!--       Calculate surface temperature
-          t_surface_h%var_1d(:)    = pt_surface * exn
-          surf_lsm_h%pt_surface(:) = pt_surface
-
 !
 !--       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
           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   
+                DO  k = nzb_soil, nzt_soil  
+                   t_soil_h%var_2d(k,m)   = soil_temperature(k)
+                   m_soil_h%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_h%m_wilt(k,m))
+                   soil_moisture(k)       = MAX(soil_moisture(k),wilting_point)           
+                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_layer(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_layer(nzb_soil)     &
+                                                 * 0.25_wp
+                surf_lsm_h%lambda_h_def(m)    = pavement_heat_conduct                                     
+                surf_lsm_h%rho_c_def(m)       = pavement_heat_capacity                                                
+                
+             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)   = MAX(soil_moisture(k),surf_lsm_h%m_wilt(k,m))
+                   soil_moisture(k)       = MAX(soil_moisture(k),wilting_point)         
+                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
+
+             ENDIF
+ 
              i   = surf_lsm_h%i(m)            
              j   = surf_lsm_h%j(m)
              k   = surf_lsm_h%k(m)
-
+!
+!--          Calculate surface temperature. In case of bare soil, the surface 
+!--          temperature must be reset to the soil temperature in the first soil 
+!--          layer
+             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
+             ELSE
+                t_surface_h%var_1d(:)    = pt_surface * exn
+                surf_lsm_h%pt_surface(:) = pt_surface  
+             ENDIF
+             
              IF ( cloud_physics )  THEN
                 pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
@@ -2079,34 +2374,82 @@
              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)
-          ENDDO
-!
-!--       Vertical surfaces
-          DO  l = 0, 3
-
-             t_soil_v(l)%var_2d    = 0.0_wp
-             m_soil_v(l)%var_2d    = 0.0_wp
-             m_liq_eb_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)
-             DO  k = nzb_soil, nzt_soil
-                t_soil_v(l)%var_2d(k,:)      = soil_temperature(k)
-                m_soil_v(l)%var_2d(k,:)      = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(:))
-                soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
-             ENDDO
-             t_soil_v(l)%var_2d(nzt_soil+1,:) = soil_temperature(nzt_soil+1)
-
-!
-!--          Calculate surface temperature
-             t_surface_v(l)%var_1d(:)   = pt_surface * exn
-             surf_lsm_h%pt_surface(:)   = pt_surface
-
-!
-!--          Set artifical values for ts and us so that r_a has its initial value 
-!--          for the first time step. Only for interior core domain, not for ghost points
+             surf_lsm_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)   = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(k,m))
+                      soil_moisture(k)          = MAX(soil_moisture(k),wilting_point)   
+                   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_layer(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_layer(nzb_soil)      &
+                                                * 0.25_wp
+                   surf_lsm_v(l)%lambda_h_def(m)   = pavement_heat_conduct                                     
+                   surf_lsm_v(l)%rho_c_def(m)      = pavement_heat_capacity                           
+                   
+                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)   = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(k,m))
+                      soil_moisture(k)          = MAX(soil_moisture(k),wilting_point)         
+                   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        
+         
+!
+!--             Calculate surface temperature. In case of bare soil, the surface 
+!--             temperature must be reset to the soil temperature in the first soil 
+!--             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
+                ELSE
+                   t_surface_v(l)%var_1d(:)      = pt_surface * exn
+                   surf_lsm_v(l)%pt_surface(:)   = pt_surface               
+                ENDIF
+
+!
+!--             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)
@@ -2121,9 +2464,9 @@
 !
 !--             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 ( 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_h%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)
+                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
              ENDDO
           ENDDO
@@ -2191,93 +2534,29 @@
        ENDDO
 
-       IF ( veg_type /= 0 )  THEN
-          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
-             min_canopy_resistance = veg_pars(0,veg_type)
-          ENDIF
-          IF ( leaf_area_index == 9999999.9_wp )  THEN
-             leaf_area_index = veg_pars(1,veg_type)         
-          ENDIF
-          IF ( vegetation_coverage == 9999999.9_wp )  THEN
-             vegetation_coverage = veg_pars(2,veg_type)     
-          ENDIF
-          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
-              canopy_resistance_coefficient= veg_pars(3,veg_type)     
-          ENDIF
-          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-             lambda_surface_stable = surface_pars(0,veg_type)         
-          ENDIF
-          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-             lambda_surface_unstable = surface_pars(1,veg_type)        
-          ENDIF
-          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
-             f_shortwave_incoming = surface_pars(2,veg_type)        
-          ENDIF
-          IF ( z0_eb == 9999999.9_wp )  THEN
-             roughness_length = roughness_par(0,veg_type) 
-             z0_eb            = roughness_par(0,veg_type) 
-          ENDIF
-          IF ( z0h_eb == 9999999.9_wp )  THEN
-             z0h_eb = roughness_par(1,veg_type)
-          ENDIF
-          IF ( z0q_eb == 9999999.9_wp )  THEN
-             z0q_eb = roughness_par(2,veg_type)
-          ENDIF
-          z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp )
-
-          IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
-             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_distribution(k,veg_type)
-                   root_fraction(k)        = root_distribution(k,veg_type)
+
+
+       IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
+          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_distribution(k,vegetation_type)
+                root_fraction(k)        = root_distribution(k,vegetation_type)
+             ENDDO
+          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
+                   surf_lsm_v(l)%root_fr(k,m) = root_distribution(k,vegetation_type)
+                   root_fraction(k)           = root_distribution(k,vegetation_type)
                 ENDDO
              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
-                      surf_lsm_v(l)%root_fr(k,m) = root_distribution(k,veg_type)
-                      root_fraction(k)    = root_distribution(k,veg_type)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ELSE
-
-          IF ( z0_eb == 9999999.9_wp )  THEN
-             z0_eb = roughness_length
-          ENDIF
-          IF ( z0h_eb == 9999999.9_wp )  THEN
-             z0h_eb = z0_eb * z0h_factor
-          ENDIF
-          IF ( z0q_eb == 9999999.9_wp )  THEN
-             z0q_eb = z0_eb * z0h_factor
-          ENDIF
-
-       ENDIF
-
-!
-!--    For surfaces covered with pavement, set depth of the pavement (with dry
-!--    soil below). The depth must be greater than the first soil layer depth
-       IF ( veg_type == 20 )  THEN
-          IF ( pave_depth == 9999999.9_wp )  THEN
-             pave_depth = zs(nzb_soil)  
-          ELSE
-             pave_depth = MAX( zs(nzb_soil), pave_depth )
-          ENDIF
-       ENDIF
-
-!
-!--    Map vegetation and soil types to 2D array to allow for heterogeneous 
-!--    surfaces via user interface see below
-!--    First for horizontal surfaces
-       surf_lsm_h%veg_type_2d  = veg_type
-       surf_lsm_h%soil_type_2d = soil_type
+           ENDDO
+       ENDIF
 
 !
@@ -2287,15 +2566,8 @@
        surf_lsm_h%c_veg             = vegetation_coverage
        surf_lsm_h%g_d               = canopy_resistance_coefficient
-       surf_lsm_h%lambda_surface_s  = lambda_surface_stable
-       surf_lsm_h%lambda_surface_u  = lambda_surface_unstable
        surf_lsm_h%f_sw_in           = f_shortwave_incoming
-       surf_lsm_h%z0                = z0_eb
-       surf_lsm_h%z0h               = z0h_eb
-       surf_lsm_h%z0q               = z0q_eb
 
 !--    Vertical surfaces
        DO  l = 0, 3
-          surf_lsm_v(l)%veg_type_2d  = veg_type
-          surf_lsm_v(l)%soil_type_2d = soil_type
 
 !
@@ -2305,34 +2577,11 @@
           surf_lsm_v(l)%c_veg             = vegetation_coverage
           surf_lsm_v(l)%g_d               = canopy_resistance_coefficient
-          surf_lsm_v(l)%lambda_surface_s  = lambda_surface_stable
-          surf_lsm_v(l)%lambda_surface_u  = lambda_surface_unstable
           surf_lsm_v(l)%f_sw_in           = f_shortwave_incoming
-          surf_lsm_v(l)%z0                = z0_eb
-          surf_lsm_v(l)%z0h               = z0h_eb
-          surf_lsm_v(l)%z0q               = z0q_eb
        ENDDO
-
+ 
 !
 !--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
        CALL user_init_land_surface
 
-!
-!--    Set flag parameter if vegetation type was set to a water surface. Also
-!--    set temperature to a constant value in all "soil" layers.
-!--    For now, do not set water surface for vertical surfaces
-       DO  m = 1, surf_lsm_h%ns
-          i   = surf_lsm_h%i(m)            
-          j   = surf_lsm_h%j(m)
-
-          IF ( surf_lsm_h%veg_type_2d(m) == 14  .OR.                           &
-               surf_lsm_h%veg_type_2d(m) == 15 )  THEN
-             surf_lsm_h%water_surface(m) = .TRUE.
-             t_soil_h%var_2d(:,m) = t_surface_h%var_1d(m)
-          ELSEIF ( surf_lsm_h%veg_type_2d(m) == 20 )  THEN
-             surf_lsm_h%pave_surface(m) = .TRUE.
-             m_soil_h%var_2d(:,m)            = 0.0_wp
-          ENDIF
-
-       ENDDO
 
 !
@@ -2343,10 +2592,10 @@
        t_soil_h_p    = t_soil_h
        m_soil_h_p    = m_soil_h
-       m_liq_eb_h_p  = m_liq_eb_h
+       m_liq_h_p  = m_liq_h
        t_surface_h_p = t_surface_h
 
        t_soil_v_p    = t_soil_v
        m_soil_v_p    = m_soil_v
-       m_liq_eb_v_p  = m_liq_eb_v
+       m_liq_v_p  = m_liq_v
        t_surface_v_p = t_surface_v
 
@@ -2377,8 +2626,8 @@
 !
 !--    Allocate global 1D arrays
-       ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
-       ALLOCATE ( ddz_soil_stag(nzb_soil:nzt_soil) )
-       ALLOCATE ( dz_soil(nzb_soil:nzt_soil+1) )
-       ALLOCATE ( dz_soil_stag(nzb_soil:nzt_soil) )
+       ALLOCATE ( ddz_soil(nzb_soil:nzt_soil) )
+       ALLOCATE ( ddz_soil_layer(nzb_soil:nzt_soil+1) )
+       ALLOCATE ( dz_soil(nzb_soil:nzt_soil) )
+       ALLOCATE ( dz_soil_layer(nzb_soil:nzt_soil+1) )
        ALLOCATE ( root_extr(nzb_soil:nzt_soil) )
   
@@ -2393,6 +2642,6 @@
 !
 !--    Horizontal surfaces
-       ALLOCATE ( m_liq_eb_h%var_1d(1:surf_lsm_h%ns)                     )
-       ALLOCATE ( m_liq_eb_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)                  )
@@ -2401,9 +2650,10 @@
        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) )
+
 !
 !--    Vertical surfaces
        DO  l = 0, 3
-          ALLOCATE ( m_liq_eb_v(l)%var_1d(1:surf_lsm_v(l)%ns)                     )
-          ALLOCATE ( m_liq_eb_v(l)_p%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(l)_p%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(l)_p%var_1d(1:surf_lsm_v(l)%ns)                  )
@@ -2416,6 +2666,6 @@
 !
 !--    Horizontal surfaces
-       ALLOCATE ( m_liq_eb_h_1%var_1d(1:surf_lsm_h%ns)                   )
-       ALLOCATE ( m_liq_eb_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)                  )
@@ -2427,6 +2677,6 @@
 !--    Vertical surfaces
        DO  l = 0, 3
-          ALLOCATE ( m_liq_eb_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
-          ALLOCATE ( m_liq_eb_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)                  )
@@ -2441,5 +2691,5 @@
 !--    Allocate intermediate timestep arrays
 !--    Horizontal surfaces
-       ALLOCATE ( tm_liq_eb_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)  )
@@ -2448,5 +2698,5 @@
 !--    Horizontal surfaces
        DO  l = 0, 3
-          ALLOCATE ( tm_liq_eb_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)  )
@@ -2463,27 +2713,22 @@
 !--    Allocate 2D vegetation model arrays
 !--    Horizontal surfaces
-       ALLOCATE ( surf_lsm_h%alpha_vg(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_eb(1:surf_lsm_h%ns)           )
-       ALLOCATE ( surf_lsm_h%gamma_w_sat(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%l_vg(1:surf_lsm_h%ns)             )
+       ALLOCATE ( surf_lsm_h%lambda_h_def(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%m_fc(1:surf_lsm_h%ns)             )
-       ALLOCATE ( surf_lsm_h%m_res(1:surf_lsm_h%ns)            )
-       ALLOCATE ( surf_lsm_h%m_sat(1:surf_lsm_h%ns)            )
-       ALLOCATE ( surf_lsm_h%m_wilt(1:surf_lsm_h%ns)           )
-       ALLOCATE ( surf_lsm_h%n_vg(1:surf_lsm_h%ns)             )
-       ALLOCATE ( surf_lsm_h%pave_surface(1:surf_lsm_h%ns)     )
-       ALLOCATE ( surf_lsm_h%qsws_eb(1:surf_lsm_h%ns)          )
-       ALLOCATE ( surf_lsm_h%qsws_soil_eb(1:surf_lsm_h%ns)     )
-       ALLOCATE ( surf_lsm_h%qsws_liq_eb(1:surf_lsm_h%ns)      )
-       ALLOCATE ( surf_lsm_h%qsws_veg_eb(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%rho_c_def(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)         )
@@ -2492,37 +2737,30 @@
        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%shf_eb(1:surf_lsm_h%ns)           )
-       ALLOCATE ( surf_lsm_h%soil_type_2d(1:surf_lsm_h%ns)     )
-       ALLOCATE ( surf_lsm_h%veg_type_2d(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%pave_surface  = .FALSE.
+       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)%alpha_vg(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_eb(1:surf_lsm_v(l)%ns)           )
-          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(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)%l_vg(1:surf_lsm_v(l)%ns)             )
+          ALLOCATE ( surf_lsm_v(l)%lambda_h_def(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)%m_fc(1:surf_lsm_v(l)%ns)             )
-          ALLOCATE ( surf_lsm_v(l)%m_res(1:surf_lsm_v(l)%ns)            )
-          ALLOCATE ( surf_lsm_v(l)%m_sat(1:surf_lsm_v(l)%ns)            )
-          ALLOCATE ( surf_lsm_v(l)%m_wilt(1:surf_lsm_v(l)%ns)           )
-          ALLOCATE ( surf_lsm_v(l)%n_vg(1:surf_lsm_v(l)%ns)             )
-          ALLOCATE ( surf_lsm_v(l)%pave_surface(1:surf_lsm_v(l)%ns)     )
-          ALLOCATE ( surf_lsm_v(l)%qsws_eb(1:surf_lsm_v(l)%ns)          )
-          ALLOCATE ( surf_lsm_v(l)%qsws_soil_eb(1:surf_lsm_v(l)%ns)     )
-          ALLOCATE ( surf_lsm_v(l)%qsws_liq_eb(1:surf_lsm_v(l)%ns)      )
-          ALLOCATE ( surf_lsm_v(l)%qsws_veg_eb(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)%rho_c_def(1:surf_lsm_h%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)         )
@@ -2531,11 +2769,10 @@
           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)%shf_eb(1:surf_lsm_v(l)%ns)           )
-          ALLOCATE ( surf_lsm_v(l)%soil_type_2d(1:surf_lsm_v(l)%ns)     )
-          ALLOCATE ( surf_lsm_v(l)%veg_type_2d(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.
-          surf_lsm_v(l)%pave_surface  = .FALSE.
+          surf_lsm_v(l)%water_surface     = .FALSE.
+          surf_lsm_v(l)%pavement_surface  = .FALSE.
+          surf_lsm_v(l)%vegetation_surface   = .FALSE. 
+          
        ENDDO
    
@@ -2547,5 +2784,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_eb_h  => m_liq_eb_h_1;  m_liq_eb_h_p  => m_liq_eb_h_2
+       m_liq_h  => m_liq_h_1;  m_liq_h_p  => m_liq_h_2
 !
 !--    Vertical surfaces
@@ -2553,5 +2790,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_eb_v  => m_liq_eb_v_1;  m_liq_eb_v_p  => m_liq_eb_v_2
+       m_liq_v  => m_liq_v_1;  m_liq_v_p  => m_liq_v_2
 #endif
 
@@ -2574,4 +2811,5 @@
        NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
                                   canopy_resistance_coefficient,               &
+                                  constant_roughness,                          &
                                   conserve_water_content,                      &
                                   f_shortwave_incoming, field_capacity,        & 
@@ -2581,11 +2819,16 @@
                                   l_vangenuchten, min_canopy_resistance,       &
                                   min_soil_resistance, n_vangenuchten,         &
-                                  pave_depth, pave_heat_capacity,              &
-                                  pave_heat_conductivity,                      &
+                                  pavement_depth, pavement_heat_capacity,      &
+                                  pavement_heat_conduct, pavement_type,        &
                                   residual_moisture, root_fraction,            &
                                   saturation_moisture, skip_time_do_lsm,       &
                                   soil_moisture, soil_temperature, soil_type,  &
-                                  vegetation_coverage, veg_type, wilting_point,& 
-                                  zs, z0_eb, z0h_eb, z0q_eb
+                                  surface_type,                                &
+                                  vegetation_coverage, vegetation_type,        &
+                                  water_temperature, water_type,               &
+                                  wilting_point, zs, z0_vegetation,            &
+                                  z0h_vegetation, z0q_vegetation, z0_water,    &
+                                  z0h_water, z0q_water, z0_pavement,           &
+                                  z0h_pavement, z0q_pavement
        
        line = ' '
@@ -2668,16 +2911,11 @@
        DO  m = 1, surf%ns
 
-          IF ( surf%pave_surface(m) )  THEN
-             surf%rho_c_total(nzb_soil,m) = pave_heat_capacity
-             lambda_temp(nzb_soil)          = pave_heat_conductivity
-          ENDIF
-
           IF (  .NOT.  surf%water_surface(m) )  THEN
              DO  k = nzb_soil, nzt_soil
 
-                IF ( surf%pave_surface(m)  .AND.  zs(k) <= pave_depth )  THEN
+                IF ( surf%pavement_surface(m)  .AND.  zs(k) <= pavement_depth )  THEN
                     
-                   surf%rho_c_total(k,m) = pave_heat_capacity
-                   lambda_temp(k)        = pave_heat_conductivity   
+                   surf%rho_c_total(k,m) = surf%rho_c_def(m)
+                   lambda_temp(k)   = surf%lambda_h_def(m) 
 
                 ELSE            
@@ -2686,5 +2924,5 @@
 !--                into account water content 
                    surf%rho_c_total(k,m) = (rho_c_soil *                 &
-                                               ( 1.0_wp - surf%m_sat(m) )&
+                                               ( 1.0_wp - surf%m_sat(k,m) )&
                                                + rho_c_water * surf_m_soil%var_2d(k,m) )
 
@@ -2692,9 +2930,9 @@
 !--                Calculate soil heat conductivity at the center of the soil
 !--                layers
-                   lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(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)
 
-                   ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(k,m) /             &
-                                                     surf%m_sat(m) ) )
+                   ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(k,m) / &
+                                                     surf%m_sat(k,m) ) )
 
                    lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +       &
@@ -2704,18 +2942,20 @@
 
 !
-!--          Calculate soil heat conductivity (lambda_h) at the _stag level 
+!--          Calculate soil heat conductivity (lambda_h) at the _layer level 
 !--          using linear interpolation. For pavement surface, the
 !--          true pavement depth is considered
              DO  k = nzb_soil, nzt_soil-1
-                IF ( surf%pave_surface(m)  .AND.  zs(k)   < pave_depth   &
-                                                 .AND.  zs(k+1) > pave_depth )  THEN
-                   surf%lambda_h(k,m) = ( pave_depth - zs(k) ) / dz_soil(k+1)&
-                                     * lambda_temp(k)                          &
-                                     + ( 1.0_wp - ( pave_depth - zs(k) )       &
-                                     / dz_soil(k+1) ) * lambda_temp(k+1)
+                IF ( surf%pavement_surface(m)  .AND.  zs(k) < pavement_depth   &
+                     .AND.  zs(k+1) > pavement_depth )                         &
+                THEN
+                   surf%lambda_h(k,m) = ( pavement_depth - zs(k) ) * ddz_soil(k+1)   &
+                                      * lambda_temp(k)                               &
+                                      + ( zs(k+1) - pavement_depth ) * ddz_soil(k+1) &
+                                      * lambda_temp(k+1)
                 ELSE
-                   surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )&
-                                     * 0.5_wp
+                   surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
+                                        * 0.5_wp
                 ENDIF
+
              ENDDO
              surf%lambda_h(nzt_soil,m) = lambda_temp(nzt_soil)
@@ -2725,18 +2965,18 @@
              tend(:) = 0.0_wp
 
-             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *&
+             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(nzb_soil+1)            &
-                      + surf%ghf_eb(m) ) * ddz_soil_stag(nzb_soil)
+                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil(nzb_soil)        &
+                      + surf%ghf(m) ) * ddz_soil_layer(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(k+1)                                  &
-                              - surf%lambda_h(k-1,m)                     &
-                              * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) )                &
-                              * ddz_soil(k)                                    &
-                            ) * ddz_soil_stag(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(k)                                             &
+                              - surf%lambda_h(k-1,m)                                    &
+                              * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
+                              * ddz_soil(k-1)                                           &
+                            ) * ddz_soil_layer(k)
 
              ENDDO
@@ -2769,40 +3009,31 @@
 !
 !--             Calculate soil diffusivity at the center of the soil layers
-                lambda_temp(k) = (- b_ch * surf%gamma_w_sat(m) * psi_sat     &
-                                  / surf%m_sat(m) ) * ( MAX( surf_m_soil%var_2d(k,m),    &
-                                  surf%m_wilt(m) ) / surf%m_sat(m) )**(&
+                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 )
 
 !
 !--             Parametrization of Van Genuchten
-                IF ( soil_type /= 7 )  THEN
-!
-!--                Calculate the hydraulic conductivity after Van Genuchten 
-!--                (1980)
-                   h_vg = ( ( ( surf%m_res(m) - surf%m_sat(m) ) /  &
-                              ( surf%m_res(m) -                          &
-                                MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(m) )       &
-                              )                                                &
-                            )**(                                               &
-                          surf%n_vg(m) / ( surf%n_vg(m) - 1.0_wp ) &
-                               ) - 1.0_wp                                      &
-                          )**( 1.0_wp / surf%n_vg(m) ) / surf%alpha_vg(m)
-
-                   gamma_temp(k) = surf%gamma_w_sat(m) * ( ( ( 1.0_wp +  &
-                          ( surf%alpha_vg(m) * h_vg )**surf%n_vg(m)&
+!--             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                                      &
+                       )**( 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(m)) - (        &
-                          surf%alpha_vg(m) * h_vg )**( surf%n_vg(m)&
-                              - 1.0_wp) )**2 )                                 &
-                              / ( ( 1.0_wp + ( surf%alpha_vg(m) * h_vg   &
-                              )**surf%n_vg(m) )**( ( 1.0_wp  - 1.0_wp    &
-                              / surf%n_vg(m) ) *                         &
-                              ( surf%l_vg(m) + 2.0_wp) ) )
-!
-!--             Parametrization of Clapp & Hornberger
-                ELSE
-                   gamma_temp(k) = surf%gamma_w_sat(m) * ( surf_m_soil%var_2d(k,m)   &
-                                   / surf%m_sat(m) )**(2.0_wp * b_ch + 3.0_wp)
-                ENDIF
+                           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
@@ -2812,7 +3043,7 @@
 !--          when humidity is enabled in the atmosphere and the surface type
 !--          is not pavement (implies dry soil below).
-             IF ( humidity  .AND.  .NOT.  surf%pave_surface(m) )  THEN
-!
-!--             Calculate soil diffusivity (lambda_w) at the _stag level 
+             IF ( humidity  .AND.  .NOT.  surf%pavement_surface(m) )  THEN
+!
+!--             Calculate soil diffusivity (lambda_w) at the _layer level 
 !--             using linear interpolation. To do: replace this with
 !--             ECMWF-IFS Eq. 8.81
@@ -2837,5 +3068,5 @@
                 ENDIF     
 
-!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l     
+!--             The root extraction (= root_extr * qsws_veg / (rho_l     
 !--             * l_v)) ensures the mass conservation for water. The         
 !--             transpiration of plants equals the cumulative withdrawals by 
@@ -2852,5 +3083,5 @@
                 m_total = 0.0_wp
                 DO  k = nzb_soil, nzt_soil
-                    IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(m) )  THEN
+                    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)
                     ENDIF
@@ -2858,5 +3089,5 @@
                  IF ( m_total > 0.0_wp )  THEN
                    DO  k = nzb_soil, nzt_soil
-                      IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(m) )  THEN
+                      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
@@ -2872,19 +3103,19 @@
                 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(nzb_soil+1) - surf%gamma_w(nzb_soil,m) - &
+                         * ddz_soil(nzb_soil) - surf%gamma_w(nzb_soil,m) - &
                                                      (                         &
-                            root_extr(nzb_soil) * surf%qsws_veg_eb(m)    &
-                            + surf%qsws_soil_eb(m) ) * drho_l_lv )       &
-                            * ddz_soil_stag(nzb_soil)
+                            root_extr(nzb_soil) * surf%qsws_veg(m)    &
+                            + surf%qsws_soil(m) ) * drho_l_lv )       &
+                            * ddz_soil_layer(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(k+1)                   &
+                             - surf_m_soil%var_2d(k,m) ) * ddz_soil(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(k)                      &
+                             surf_m_soil%var_2d(k-1,m)) * ddz_soil(k-1)                      &
                              + surf%gamma_w(k-1,m) - (root_extr(k)       &
-                             * surf%qsws_veg_eb(m) * drho_l_lv)          &
-                             ) * ddz_soil_stag(k)
+                             * surf%qsws_veg(m) * drho_l_lv)          &
+                             ) * ddz_soil_layer(k)
                 ENDDO
                 tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)            &
@@ -2892,9 +3123,9 @@
                                      * ( surf_m_soil%var_2d(nzt_soil,m)                    &
                                      - surf_m_soil%var_2d(nzt_soil-1,m))                   &
-                                     * ddz_soil(nzt_soil)                      &
+                                     * ddz_soil(nzt_soil-1)                      &
                                      + surf%gamma_w(nzt_soil-1,m) - (    &
                                        root_extr(nzt_soil)                     &
-                                     * surf%qsws_veg_eb(m) * drho_l_lv ) &
-                                  ) * ddz_soil_stag(nzt_soil)             
+                                     * surf%qsws_veg(m) * drho_l_lv ) &
+                                  ) * ddz_soil_layer(nzt_soil)             
 
                 surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_m_soil%var_2d(nzb_soil:nzt_soil,m)    &
@@ -2950,5 +3181,5 @@
        IF ( humidity )  THEN
           m_soil_h    = m_soil_h_p
-          m_liq_eb_h  = m_liq_eb_h_p
+          m_liq_h  = m_liq_h_p
        ENDIF
 !
@@ -2958,5 +3189,5 @@
        IF ( humidity )  THEN
           m_soil_v    = m_soil_v_p
-          m_liq_eb_v  = m_liq_eb_v_p
+          m_liq_v  = m_liq_v_p
        ENDIF
 
@@ -2972,5 +3203,5 @@
              IF ( humidity )  THEN
                 m_soil_h    => m_soil_h_1;   m_soil_h_p    => m_soil_h_2
-                m_liq_eb_h  => m_liq_eb_h_1; m_liq_eb_h_p  => m_liq_eb_h_2
+                m_liq_h  => m_liq_h_1; m_liq_h_p  => m_liq_h_2
              ENDIF
 !
@@ -2980,5 +3211,5 @@
              IF ( humidity )  THEN
                 m_soil_v    => m_soil_v_1;   m_soil_v_p    => m_soil_v_2
-                m_liq_eb_v  => m_liq_eb_v_1; m_liq_eb_v_p  => m_liq_eb_v_2
+                m_liq_v  => m_liq_v_1; m_liq_v_p  => m_liq_v_2
              ENDIF
 
@@ -2992,5 +3223,5 @@
              IF ( humidity )  THEN
                 m_soil_h    => m_soil_h_2;   m_soil_h_p    => m_soil_h_1
-                m_liq_eb_h  => m_liq_eb_h_2; m_liq_eb_h_p  => m_liq_eb_h_1
+                m_liq_h  => m_liq_h_2; m_liq_h_p  => m_liq_h_1
              ENDIF
 !
@@ -3000,5 +3231,5 @@
              IF ( humidity )  THEN
                 m_soil_v    => m_soil_v_2;   m_soil_v_p    => m_soil_v_1
-                m_liq_eb_v  => m_liq_eb_v_2; m_liq_eb_v_p  => m_liq_eb_v_1
+                m_liq_v  => m_liq_v_2; m_liq_v_p  => m_liq_v_1
              ENDIF
 
@@ -3058,9 +3289,9 @@
                 c_veg_av = 0.0_wp
 
-             CASE ( 'ghf_eb*' )
-                IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
-                   ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
+             CASE ( 'ghf*' )
+                IF ( .NOT. ALLOCATED( ghf_av ) )  THEN
+                   ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
                 ENDIF
-                ghf_eb_av = 0.0_wp
+                ghf_av = 0.0_wp
 
              CASE ( 'lai*' )
@@ -3070,9 +3301,9 @@
                 lai_av = 0.0_wp
 
-             CASE ( 'm_liq_eb*' )
-                IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
-                   ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+             CASE ( 'm_liq*' )
+                IF ( .NOT. ALLOCATED( m_liq_av ) )  THEN
+                   ALLOCATE( m_liq_av(nysg:nyng,nxlg:nxrg) )
                 ENDIF
-                m_liq_eb_av = 0.0_wp
+                m_liq_av = 0.0_wp
 
              CASE ( 'm_soil' )
@@ -3082,27 +3313,21 @@
                 m_soil_av = 0.0_wp
 
-             CASE ( 'qsws_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
-                   ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
+             CASE ( 'qsws_liq*' )
+                IF ( .NOT. ALLOCATED( qsws_liq_av ) )  THEN
+                   ALLOCATE( qsws_liq_av(nysg:nyng,nxlg:nxrg) )
                 ENDIF
-                qsws_eb_av = 0.0_wp
-
-             CASE ( 'qsws_liq_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
-                   ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                qsws_liq_av = 0.0_wp
+
+             CASE ( 'qsws_soil*' )
+                IF ( .NOT. ALLOCATED( qsws_soil_av ) )  THEN
+                   ALLOCATE( qsws_soil_av(nysg:nyng,nxlg:nxrg) )
                 ENDIF
-                qsws_liq_eb_av = 0.0_wp
-
-             CASE ( 'qsws_soil_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
-                   ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
+                qsws_soil_av = 0.0_wp
+
+             CASE ( 'qsws_veg*' )
+                IF ( .NOT. ALLOCATED( qsws_veg_av ) )  THEN
+                   ALLOCATE( qsws_veg_av(nysg:nyng,nxlg:nxrg) )
                 ENDIF
-                qsws_soil_eb_av = 0.0_wp
-
-             CASE ( 'qsws_veg_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
-                   ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                qsws_veg_eb_av = 0.0_wp
+                qsws_veg_av = 0.0_wp
 
              CASE ( 'r_a*' )
@@ -3118,10 +3343,4 @@
                 r_s_av = 0.0_wp
 
-             CASE ( 'shf_eb*' )
-                IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
-                   ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                shf_eb_av = 0.0_wp
-
              CASE ( 't_soil' )
                 IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
@@ -3160,9 +3379,9 @@
              ENDDO
 
-          CASE ( 'ghf_eb*' )
+          CASE ( 'ghf*' )
              DO  m = 1, surf_lsm_h%ns
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                ghf_eb_av(j,i) = ghf_eb_av(j,i) + surf_lsm_h%ghf_eb(m)
+                ghf_av(j,i) = ghf_av(j,i) + surf_lsm_h%ghf(m)
              ENDDO
 
@@ -3174,9 +3393,9 @@
              ENDDO
 
-          CASE ( 'm_liq_eb*' )
+          CASE ( 'm_liq*' )
              DO  m = 1, surf_lsm_h%ns
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb_h%var_1d(m)
+                m_liq_av(j,i) = m_liq_av(j,i) + m_liq_h%var_1d(m)
              ENDDO
 
@@ -3190,33 +3409,26 @@
              ENDDO
 
-          CASE ( 'qsws_eb*' )
+          CASE ( 'qsws_liq*' )
              DO  m = 1, surf_lsm_h%ns
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                qsws_eb_av(j,i) = qsws_eb_av(j,i) + surf_lsm_h%qsws_eb(m)
-             ENDDO
-
-          CASE ( 'qsws_liq_eb*' )
+                qsws_liq_av(j,i) = qsws_liq_av(j,i) +                    &
+                                      surf_lsm_h%qsws_liq(m)
+             ENDDO
+
+          CASE ( 'qsws_soil*' )
              DO  m = 1, surf_lsm_h%ns
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) +                    &
-                                      surf_lsm_h%qsws_liq_eb(m)
-             ENDDO
-
-          CASE ( 'qsws_soil_eb*' )
+                qsws_soil_av(j,i) = qsws_soil_av(j,i) +                  &
+                                       surf_lsm_h%qsws_soil(m)
+             ENDDO
+
+          CASE ( 'qsws_veg*' )
              DO  m = 1, surf_lsm_h%ns
                 i   = surf_lsm_h%i(m)            
                 j   = surf_lsm_h%j(m)
-                qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) +                  &
-                                       surf_lsm_h%qsws_soil_eb(m)
-             ENDDO
-
-          CASE ( 'qsws_veg_eb*' )
-             DO  m = 1, surf_lsm_h%ns
-                i   = surf_lsm_h%i(m)            
-                j   = surf_lsm_h%j(m)
-                qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) +                    &
-                                      surf_lsm_h%qsws_veg_eb(m)
+                qsws_veg_av(j,i) = qsws_veg_av(j,i) +                    &
+                                      surf_lsm_h%qsws_veg(m)
              ENDDO
 
@@ -3233,11 +3445,4 @@
                 j   = surf_lsm_h%j(m)
                 r_s_av(j,i) = r_s_av(j,i) + surf_lsm_h%r_s(m)
-             ENDDO
-
-          CASE ( 'shf_eb*' )
-             DO  m = 1, surf_lsm_h%ns
-                i   = surf_lsm_h%i(m)            
-                j   = surf_lsm_h%j(m)
-                shf_eb_av(j,i) = shf_eb_av(j,i) + surf_lsm_h%shf_eb(m)
              ENDDO
 
@@ -3281,8 +3486,8 @@
              ENDDO
 
-          CASE ( 'ghf_eb*' )
+          CASE ( 'ghf*' )
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   ghf_eb_av(j,i) = ghf_eb_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
@@ -3295,8 +3500,8 @@
              ENDDO
 
-          CASE ( 'm_liq_eb*' )
+          CASE ( 'm_liq*' )
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   m_liq_eb_av(j,i) = m_liq_eb_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
@@ -3311,29 +3516,22 @@
              ENDDO
 
-          CASE ( 'qsws_eb*' )
+          CASE ( 'qsws_liq*' )
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   qsws_eb_av(j,i) = qsws_eb_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
 
-          CASE ( 'qsws_liq_eb*' )
+          CASE ( 'qsws_soil*' )
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   qsws_liq_eb_av(j,i) = qsws_liq_eb_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
 
-          CASE ( 'qsws_soil_eb*' )
+          CASE ( 'qsws_veg*' )
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-          CASE ( 'qsws_veg_eb*' )
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   qsws_veg_eb_av(j,i) = qsws_veg_eb_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
@@ -3503,15 +3701,15 @@
           grid = 'zu1'
 
-       CASE ( 'ghf_eb*_xy' )        ! 2d-array
+       CASE ( 'ghf*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
                 i                   = surf_lsm_h%i(m)            
                 j                   = surf_lsm_h%j(m)
-                local_pf(i,j,nzb+1) = surf_lsm_h%ghf_eb(m)
+                local_pf(i,j,nzb+1) = surf_lsm_h%ghf(m)
              ENDDO
           ELSE
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
+                   local_pf(i,j,nzb+1) = ghf_av(j,i)
                 ENDDO
              ENDDO
@@ -3539,15 +3737,15 @@
           grid = 'zu1'
 
-       CASE ( 'm_liq_eb*_xy' )        ! 2d-array
+       CASE ( 'm_liq*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
              DO  m = 1, surf_lsm_h%ns
                 i                   = surf_lsm_h%i(m)            
                 j                   = surf_lsm_h%j(m)
-                local_pf(i,j,nzb+1) = m_liq_eb_h%var_1d(m)
+                local_pf(i,j,nzb+1) = m_liq_h%var_1d(m)
              ENDDO
           ELSE
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
+                   local_pf(i,j,nzb+1) = m_liq_av(j,i)
                 ENDDO
              ENDDO
@@ -3581,15 +3779,15 @@
           IF ( mode == 'xy' ) grid = 'zs'
 
-       CASE ( 'qsws_eb*_xy' )        ! 2d-array
+       CASE ( 'qsws_liq*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
              DO  m = 1, surf_lsm_h%ns
                 i                   = surf_lsm_h%i(m)            
                 j                   = surf_lsm_h%j(m)
-                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_eb(m)
+                local_pf(i,j,nzb+1) = surf_lsm_h%qsws_liq(m)
              ENDDO
           ELSE
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng 
-                   local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
+                   local_pf(i,j,nzb+1) =  qsws_liq_av(j,i)
                 ENDDO
              ENDDO
@@ -3599,15 +3797,15 @@
           grid = 'zu1'
 
-       CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
+       CASE ( 'qsws_soil*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
              DO  m = 1, surf_lsm_h%ns
                 i                   = surf_lsm_h%i(m)            
                 j                   = surf_lsm_h%j(m)
-                local_pf(i,j,nzb+1) = surf_lsm_h%qsws_liq_eb(m)
+                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_soil(m)
              ENDDO
           ELSE
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng 
-                   local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
+                   local_pf(i,j,nzb+1) =  qsws_soil_av(j,i)
                 ENDDO
              ENDDO
@@ -3617,33 +3815,15 @@
           grid = 'zu1'
 
-       CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
+       CASE ( 'qsws_veg*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
              DO  m = 1, surf_lsm_h%ns
                 i                   = surf_lsm_h%i(m)            
                 j                   = surf_lsm_h%j(m)
-                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_soil_eb(m)
+                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_veg(m)
              ENDDO
           ELSE
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng 
-                   local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
-                ENDDO
-             ENDDO
-          ENDIF
-
-          two_d = .TRUE.
-          grid = 'zu1'
-
-       CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
-          IF ( av == 0 ) THEN
-             DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
-                j                   = surf_lsm_h%j(m)
-                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_veg_eb(m)
-             ENDDO
-          ELSE
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng 
-                   local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
+                   local_pf(i,j,nzb+1) =  qsws_veg_av(j,i)
                 ENDDO
              ENDDO
@@ -3683,22 +3863,4 @@
                 DO  j = nysg, nyng
                    local_pf(i,j,nzb+1) = r_s_av(j,i)
-                ENDDO
-             ENDDO
-          ENDIF
-
-          two_d = .TRUE.
-          grid = 'zu1'
-
-       CASE ( 'shf_eb*_xy' )        ! 2d-array
-          IF ( av == 0 ) THEN
-             DO  m = 1, surf_lsm_h%ns
-                i                   = surf_lsm_h%i(m)            
-                j                   = surf_lsm_h%j(m)
-                local_pf(i,j,nzb+1) =  surf_lsm_h%shf_eb(m)
-             ENDDO
-          ELSE
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
                 ENDDO
              ENDDO
@@ -3852,13 +4014,13 @@
           WRITE ( 14 )  'c_veg_av            ';  WRITE ( 14 ) c_veg_av
        ENDIF
-       IF ( ALLOCATED( ghf_eb_av ) )  THEN
-          WRITE ( 14 )  'ghf_eb_av           ';  WRITE ( 14 )  ghf_eb_av
+       IF ( ALLOCATED( ghf_av ) )  THEN
+          WRITE ( 14 )  'ghf_av           ';  WRITE ( 14 )  ghf_av
        ENDIF
        IF ( ALLOCATED( lai_av ) )  THEN
           WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
        ENDIF
-       WRITE ( 14 )     'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb_h%var_1d
-       IF ( ALLOCATED( m_liq_eb_av ) )  THEN
-          WRITE ( 14 )  'm_liq_eb_av         ';  WRITE ( 14 )  m_liq_eb_av
+       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
@@ -3866,19 +4028,13 @@
           WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
        ENDIF
-       IF ( ALLOCATED( qsws_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_eb_av          ';  WRITE ( 14 )  qsws_eb_av
-       ENDIF   
-       IF ( ALLOCATED( qsws_liq_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_liq_eb_av      ';  WRITE ( 14 )  qsws_liq_eb_av
+       IF ( ALLOCATED( qsws_liq_av ) )  THEN
+          WRITE ( 14 )  'qsws_liq_av      ';  WRITE ( 14 )  qsws_liq_av
        ENDIF  
-       IF ( ALLOCATED( qsws_soil_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_soil_eb_av     ';  WRITE ( 14 )  qsws_soil_eb_av
+       IF ( ALLOCATED( qsws_soil_av ) )  THEN
+          WRITE ( 14 )  'qsws_soil_av     ';  WRITE ( 14 )  qsws_soil_av
        ENDIF 
-       IF ( ALLOCATED( qsws_veg_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_veg_eb_av      ';  WRITE ( 14 )  qsws_veg_eb_av
+       IF ( ALLOCATED( qsws_veg_av ) )  THEN
+          WRITE ( 14 )  'qsws_veg_av      ';  WRITE ( 14 )  qsws_veg_av
        ENDIF 
-       IF ( ALLOCATED( shf_eb_av ) )  THEN
-          WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
-       ENDIF
        WRITE ( 14 )     't_soil              ';  WRITE ( 14 )  t_soil_h%var_2d
        IF ( ALLOCATED( t_soil_av ) )  THEN
@@ -4002,15 +4158,15 @@
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
-                CASE ( 'ghf_eb_av' )
-                   IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
-                      ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
+                CASE ( 'ghf_av' )
+                   IF ( .NOT. ALLOCATED( ghf_av ) )  THEN
+                      ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   ghf_eb_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_eb' )
+                CASE ( 'm_liq' )
                    IF ( k == 1 )  READ ( 13 )  tmp_walltype_1d !tmp_2d
-                   m_liq_eb_h%var_1d(1:surf_lsm_h%ns)  =                       &
+                   m_liq_h%var_1d(1:surf_lsm_h%ns)  =                       &
                                  tmp_walltype_1d(1:surf_lsm_h%ns)
 
@@ -4023,10 +4179,10 @@
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
-                CASE ( 'm_liq_eb_av' )
-                   IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
-                      ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                CASE ( 'm_liq_av' )
+                   IF ( .NOT. ALLOCATED( m_liq_av ) )  THEN
+                      ALLOCATE( m_liq_av(nysg:nyng,nxlg:nxrg) )
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   m_liq_eb_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)
 
@@ -4045,42 +4201,26 @@
                                     -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
-                CASE ( 'qsws_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
-                      ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
+                CASE ( 'qsws_liq_av' )
+                   IF ( .NOT. ALLOCATED( qsws_liq_av ) )  THEN
+                      ALLOCATE( qsws_liq_av(nysg:nyng,nxlg:nxrg) )
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_eb_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_liq_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
-                      ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                CASE ( 'qsws_soil_av' )
+                   IF ( .NOT. ALLOCATED( qsws_soil_av ) )  THEN
+                      ALLOCATE( qsws_soil_av(nysg:nyng,nxlg:nxrg) )
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_liq_eb_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)
-                CASE ( 'qsws_soil_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
-                      ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
+
+                CASE ( 'qsws_veg_av' )
+                   IF ( .NOT. ALLOCATED( qsws_veg_av ) )  THEN
+                      ALLOCATE( qsws_veg_av(nysg:nyng,nxlg:nxrg) )
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_soil_eb_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)
-
-                CASE ( 'qsws_veg_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
-                      ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF  
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_veg_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'shf_eb_av' )
-                   IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
-                      ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =       &
-                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 't_soil' )
Index: /palm/trunk/SOURCE/lpm_splitting.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_splitting.f90	(revision 2269)
+++ /palm/trunk/SOURCE/lpm_splitting.f90	(revision 2270)
@@ -25,4 +25,7 @@
 ! -----------------
 !
+! 
+! 
+! 
 ! Initial revision
 ! 
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2269)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2270)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removed 2 timeseries (shf_eb + qsws_eb). Removed _eb suffixes
+! 
+! 2265 2017-06-08 16:58:28Z schwenkel
 ! Unused variables removed.
 ! 
@@ -250,7 +253,7 @@
              'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
 
-    INTEGER(iwp) ::  dots_num  = 33  !< number of timeseries defined by default
+    INTEGER(iwp) ::  dots_num  = 25  !< number of timeseries defined by default
     INTEGER(iwp) ::  dots_soil = 26  !< starting index for soil-timeseries
-    INTEGER(iwp) ::  dots_rad  = 34  !< starting index for radiation-timeseries
+    INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries
 
     CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
@@ -263,11 +266,11 @@
              'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
              'ol           ', 'q*           ', 'w"s"         ',                &
-             's*           ', 'ghf_eb       ', 'shf_eb       ',                &
-             'qsws_eb      ', 'qsws_liq_eb  ', 'qsws_soil_eb ',                &
-             'qsws_veg_eb  ', 'r_a          ', 'r_s          ',                &
+             's*           ', 'ghf          ', 'qsws_liq     ',                &
+             'qsws_soil    ', 'qsws_veg     ', 'r_a          ',                &
+             'r_s          ',                                                  &
              'rad_net      ', 'rad_lw_in    ', 'rad_lw_out   ',                &
              'rad_sw_in    ', 'rad_sw_out   ', 'rrtm_aldif   ',                &
              'rrtm_aldir   ', 'rrtm_asdif   ', 'rrtm_asdir   ',                &                                               
-             ( 'unknown      ', i9 = 1, dots_max-42 ) /)
+             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
 
     CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
@@ -281,10 +284,10 @@
              'm            ', 'kg/kg        ', 'kg m/(kg s)  ',                &
              'kg/kg        ', 'W/m2         ', 'W/m2         ',                &
-             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
-             'W/m2         ', 's/m          ', 's/m          ',                &
+             'W/m2         ', 'W/m2         ', 's/m          ',                &
+             's/m          ',                                                  &
              'W/m2         ', 'W/m2         ', 'W/m2         ',                &
              'W/m2         ', 'W/m2         ', '             ',                &
              '             ', '             ', '             ',                &
-             ( 'unknown      ', i9 = 1, dots_max-42 ) /)
+             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
 
     CHARACTER (LEN=16) :: heatflux_output_unit     !< unit for heatflux output
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2269)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2270)
@@ -20,5 +20,4 @@
 ! Current revisions:
 ! -----------------
-! 
 ! 
 ! 
@@ -26,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Numbering changed (2 timeseries removed)
+! 
+! 2249 2017-06-06 13:58:01Z sward
 ! Allow for RRTMG runs without humidity/cloud physics
 ! 
@@ -658,12 +660,42 @@
                 CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
              ELSE
+                dopr_index(var_count) = 99
+                dopr_unit  = 'W/m2'
+                hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'rad_lw_in' )
+             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 100
+                dopr_unit  = 'W/m2'
+                hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit 
+             ENDIF
+
+          CASE ( 'rad_lw_out' )
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
                 dopr_index(var_count) = 101
                 dopr_unit  = 'W/m2'
                 hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
-                unit = dopr_unit
+                unit = dopr_unit   
              ENDIF
 
-          CASE ( 'rad_lw_in' )
-             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
+          CASE ( 'rad_sw_in' )
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
              THEN
                 message_string = 'data_output_pr = ' //                        &
@@ -676,9 +708,9 @@
                 dopr_unit  = 'W/m2'
                 hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
-                unit = dopr_unit 
+                unit = dopr_unit
              ENDIF
 
-          CASE ( 'rad_lw_out' )
-             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+          CASE ( 'rad_sw_out')
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
              THEN
                 message_string = 'data_output_pr = ' //                        &
@@ -691,38 +723,38 @@
                 dopr_unit  = 'W/m2'
                 hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
-                unit = dopr_unit   
+                unit = dopr_unit
              ENDIF
 
-          CASE ( 'rad_sw_in' )
-             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+          CASE ( 'rad_lw_cs_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
              THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(var_count) ) // ' is' // &
                                  'not available for radiation = .FALSE. or ' //&
-                                 'radiation_scheme = "constant"'
-                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
              ELSE
                 dopr_index(var_count) = 104
-                dopr_unit  = 'W/m2'
-                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                dopr_unit  = 'K/h'
+                hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
                 unit = dopr_unit
              ENDIF
 
-          CASE ( 'rad_sw_out')
-             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
+          CASE ( 'rad_lw_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
              THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(var_count) ) // ' is' // &
                                  'not available for radiation = .FALSE. or ' //&
-                                 'radiation_scheme = "constant"'
-                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
              ELSE
                 dopr_index(var_count) = 105
-                dopr_unit  = 'W/m2'
-                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                dopr_unit  = 'K/h'
+                hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
                 unit = dopr_unit
              ENDIF
 
-          CASE ( 'rad_lw_cs_hr' )
+          CASE ( 'rad_sw_cs_hr' )
              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
              THEN
@@ -739,5 +771,5 @@
              ENDIF
 
-          CASE ( 'rad_lw_hr' )
+          CASE ( 'rad_sw_hr' )
              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
              THEN
@@ -751,34 +783,4 @@
                 dopr_unit  = 'K/h'
                 hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                unit = dopr_unit
-             ENDIF
-
-          CASE ( 'rad_sw_cs_hr' )
-             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
-             THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(var_count) ) // ' is' // &
-                                 'not available for radiation = .FALSE. or ' //&
-                                 'radiation_scheme /= "rrtmg"'
-                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
-             ELSE
-                dopr_index(var_count) = 108
-                dopr_unit  = 'K/h'
-                hom(:,2,108,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                unit = dopr_unit
-             ENDIF
-
-          CASE ( 'rad_sw_hr' )
-             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
-             THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(var_count) ) // ' is' // &
-                                 'not available for radiation = .FALSE. or ' //&
-                                 'radiation_scheme /= "rrtmg"'
-                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
-             ELSE
-                dopr_index(var_count) = 109
-                dopr_unit  = 'K/h'
-                hom(:,2,109,:)  = SPREAD( zu, 2, statistic_regions+1 )
                 unit = dopr_unit
              ENDIF
Index: /palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_mod.f90	(revision 2269)
+++ /palm/trunk/SOURCE/surface_mod.f90	(revision 2270)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Parameters removed/added due to changes in the LSM
+! 
+! 2269 2017-06-09 11:57:32Z suehring
 ! Formatting and description adjustments
 ! 
@@ -137,34 +140,22 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l         !< net radiation 
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil (W/m/K) 
-
-!
-!--    Variables required for land-surface model    
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  soil_type_2d    !< soil type, 0: user-defined, 1-7: generic (see list)
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  veg_type_2d     !< vegetation type, 0: user-defined, 1-19: generic (see list)
-
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)   
+       
        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  building_surface    !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
-       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  pave_surface        !< flag parameter for pavements (asphalt etc.) (class 20)
-       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  water_surface       !< flag parameter for water surfaces (classes 14+15)
-
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  alpha_vg            !< coef. of Van Genuchten
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  pavement_surface    !< flag parameter for pavements
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  water_surface       !< flag parameter for water surfaces
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  vegetation_surface     !< flag parameter for natural land surfaces
+
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_liq               !< liquid water coverage (of vegetated area)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_veg               !< vegetation coverage
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  f_sw_in             !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ghf_eb              !< ground heat flux
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  gamma_w_sat         !< hydraulic conductivity at saturation
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ghf              !< ground heat flux
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  g_d                 !< coefficient for dependence of r_canopy on water vapour pressure deficit
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_vg                !< coef. of Van Genuchten
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  lai                 !< leaf area index
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_u    !< coupling between surface and soil (depends on vegetation type)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_s    !< coupling between surface and soil (depends on vegetation type)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  m_fc                !< soil moisture at field capacity (m3/m3)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  m_res               !< residual soil moisture
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  m_sat               !< saturation soil moisture (m3/m3)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  m_wilt              !< soil moisture at permanent wilting point (m3/m3)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  n_vg                !< coef. Van Genuchten  
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_eb             !< surface flux of latent heat (total)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_liq_eb         !< surface flux of latent heat (liquid water portion)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_soil_eb        !< surface flux of latent heat (soil portion)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg_eb         !< surface flux of latent heat (vegetation portion)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_u    !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_s    !< coupling between surface and soil (depends on vegetation type) (W/m2/K)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_liq         !< surface flux of latent heat (liquid water portion)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_soil        !< surface flux of latent heat (soil portion)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg         !< surface flux of latent heat (vegetation portion)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a                 !< aerodynamic resistance 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy            !< canopy resistance
@@ -173,9 +164,17 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_s                 !< total surface resistance (combination of r_soil and r_canopy)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy_min        !< minimum canopy (stomatal) resistance
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf_eb              !< surface flux of sensible heat
-
+
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  alpha_vg          !< coef. of Van Genuchten
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_w          !< hydraulic diffusivity of soil (?)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w           !< hydraulic conductivity of soil (W/m/K)
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total       !< volumetric heat capacity of the actual soil matrix (?) 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w_sat       !< hydraulic conductivity at saturation
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  l_vg              !< coef. of Van Genuchten
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_fc              !< soil moisture at field capacity (m3/m3)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_res             !< residual soil moisture
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_sat             !< saturation soil moisture (m3/m3)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_wilt            !< soil moisture at permanent wilting point (m3/m3)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_vg              !< coef. Van Genuchten  
+       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  rho_c_def         !< default volumetric heat capacity of the (soil) layer (J/m3/K)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total       !< volumetric heat capacity of the actual soil matrix (J/m3/K)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  root_fr           !< root fraction within the soil layers
 !
@@ -186,7 +185,7 @@
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  albedo_surf         !< albedo of the surface
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface           !< heat capacity of the wall surface skin ( J mâ2 Kâ1 )
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface           !< heat capacity of the wall surface skin (J/m2/K)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  emiss_surf          !< emissivity of the wall surface
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf         !< heat conductivity Î»S between air and surface ( W mâ2 Kâ1 )
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf         !< heat conductivity between air and surface (W/m2/K)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  roughness_wall      !< roughness relative to concrete
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_wall      !< thickness of the wall, roof and soil layers
