Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2150)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2232)
@@ -20,5 +20,10 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments to new topography and surface concept
+!   - now, also vertical walls are possible
+!   - for vertical walls, parametrization of r_a (aerodynamic resisistance) is 
+!     implemented. 
+!
+! Add check for soil moisture, it must not exceed its saturation value.
 ! 
 ! Former revisions:
@@ -159,4 +164,8 @@
 !> DALES and UCLA-LES models.
 !>
+!> @todo Extensive verification energy-balance solver for vertical surfaces, 
+!>       e.g. parametrization of r_a
+!> @todo Revise single land-surface processes for vertical surfaces, e.g. 
+!>       treatment of humidity, etc. 
 !> @todo Consider partial absorption of the net shortwave radiation by the 
 !>       skin layer.
@@ -174,6 +183,5 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  hyp, ol, pt, pt_p, prr, q, q_p, ql, qsws, shf, ts, us, vpt, z0, &
-               z0h, z0q
+        ONLY:  hyp, pt, pt_p, prr, q, q_p, ql, vpt, u, v, w
 
     USE cloud_parameters,                                                      &
@@ -183,10 +191,10 @@
         ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
                initializing_actions, intermediate_timestep_count_max,          &
-               max_masks, precipitation, pt_surface,                           &
+               land_surface, max_masks, precipitation, pt_surface,             &
                rho_surface, roughness_length, surface_pressure,                &
                timestep_scheme, tsc, z0h_factor, time_since_reference_point
 
     USE indices,                                                               &
-        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
+        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb 
 
     USE kinds
@@ -196,11 +204,18 @@
     USE radiation_model_mod,                                                   &
         ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
-               rad_lw_out_change_0, unscheduled_radiation_calls
+               rad_lw_out_change_0, radiation_scheme, unscheduled_radiation_calls
         
     USE statistics,                                                            &
         ONLY:  hom, statistic_regions
 
+    USE surface_mod,                                                        &
+        ONLY :  surf_lsm_h, surf_lsm_v, surf_type
+
     IMPLICIT NONE
 
+    TYPE surf_type_lsm
+       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1D !< 1D prognostic variable
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2D !< 2D prognostic variable
+    END TYPE surf_type_lsm
 !
 !-- LSM model constants
@@ -225,14 +240,7 @@
                     soil_type = 3    !< NAMELIST soil_type_2d
 
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  soil_type_2d, &  !< soil type, 0: user-defined, 1-7: generic (see list)
-                                                  veg_type_2d      !< vegetation type, 0: user-defined, 1-19: generic (see list)
-
-    LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface,     & !< flag parameter for water surfaces (classes 14+15)
-                                            pave_surface,      & !< flag parameter for pavements (asphalt etc.) (class 20)
-                                            building_surface     !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
-
     LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
-               force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
-               land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
+               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
 
 !   value 9999999.9_wp -> generic available or user-defined value must be set
@@ -288,100 +296,114 @@
               dz_soil                                                            !< soil grid spacing (edge-edge)
 
+
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !< surface temperature (K)
-                                                     t_surface_p, & !< progn. surface temperature (K)
-                                                     m_liq_eb,    & !< liquid water reservoir (m)
-                                                     m_liq_eb_av, & !< liquid water reservoir (m)
-                                                     m_liq_eb_p     !< progn. liquid water reservoir (m)
+    TYPE(surf_type_lsm), TARGET  ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
+                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
+                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
+                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements
+
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET  ::  &
+                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
+                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
+                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
+                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements
 #else
-    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
-                                         t_surface_p,    & 
-                                         m_liq_eb,       & 
-                                         m_liq_eb_p
-
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
-                                                     m_liq_eb_av,              &
-                                                     m_liq_eb_1, m_liq_eb_2
+    TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
+                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
+                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
+                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements 
+
+    TYPE(surf_type_lsm), TARGET  ::  t_soil_h_1,  & !<
+                                     t_soil_h_2,  & !<
+                                     m_soil_h_1,  & !<
+                                     m_soil_h_2     !<
+
+    TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
+                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
+                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
+                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
+                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements   
+
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
+                                     t_soil_v_1,  & !<
+                                     t_soil_v_2,  & !<
+                                     m_soil_v_1,  & !<
+                                     m_soil_v_2     !<
+#endif    
+
+#if defined( __nopointer )
+    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 
+
+    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 
+#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 
+
+    TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
+                                      t_surface_h_2,  & !<
+                                      m_liq_eb_h_1,   & !<
+                                      m_liq_eb_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 
+
+    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      !<
 #endif
 
-!
-!-- Temporal tendencies for time stepping            
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !< surface temperature tendency (K)
-                                             tm_liq_eb_m      !< liquid water reservoir tendency (m)
+#if defined( __nopointer )
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
+#else
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
+#endif
+
+#if defined( __nopointer )
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
+                                                        m_soil_av    !< Average of m_soil
+#else
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
+                                                        m_soil_av    !< Average of m_soil
+#endif
+
+    TYPE(surf_type_lsm), TARGET ::  tm_liq_eb_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 ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements 
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements 
+    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_soil_v_m     !< m_soil storage array, vertical surface elements 
 
 !
 !-- Energy balance variables                
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
-              alpha_vg,         & !< coef. of Van Genuchten
-              c_liq,            & !< liquid water coverage (of vegetated area)
               c_liq_av,         & !< average of c_liq
               c_soil_av,        & !< average of c_soil
-              c_veg,            & !< vegetation coverage
               c_veg_av,         & !< average of c_veg
-              f_sw_in,          & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
-              ghf_eb,           & !< ground heat flux
               ghf_eb_av,        & !< average of ghf_eb
-              gamma_w_sat,      & !< hydraulic conductivity at saturation
-              g_d,              & !< coefficient for dependence of r_canopy on water vapour pressure deficit
-              lai,              & !< leaf area index
               lai_av,           & !< average of lai
-              lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type)
-              lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type)
-              l_vg,             & !< coef. of Van Genuchten
-              m_fc,             & !< soil moisture at field capacity (m3/m3)
-              m_res,            & !< residual soil moisture
-              m_sat,            & !< saturation soil moisture (m3/m3)
-              m_wilt,           & !< soil moisture at permanent wilting point (m3/m3)
-              n_vg,             & !< coef. Van Genuchten  
-              qsws_eb,          & !< surface flux of latent heat (total)
               qsws_eb_av,       & !< average of qsws_eb
-              qsws_liq_eb,      & !< surface flux of latent heat (liquid water portion)
               qsws_liq_eb_av,   & !< average of qsws_liq_eb
-              qsws_soil_eb,     & !< surface flux of latent heat (soil portion)
               qsws_soil_eb_av,  & !< average of qsws_soil_eb
-              qsws_veg_eb,      & !< surface flux of latent heat (vegetation portion)
               qsws_veg_eb_av,   & !< average of qsws_veg_eb
-              rad_net_l,        & !< local copy of rad_net (net radiation at surface)
-              r_a,              & !< aerodynamic resistance 
               r_a_av,           & !< average of r_a
-              r_canopy,         & !< canopy resistance
-              r_soil,           & !< soil resistance
-              r_soil_min,       & !< minimum soil resistance
-              r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
               r_s_av,           & !< average of r_s
-              r_canopy_min,     & !< minimum canopy (stomatal) resistance
-              shf_eb,           & !< surface flux of sensible heat
               shf_eb_av           !< average of shf_eb
-
-
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
-              lambda_h, &   !< heat conductivity of soil (W/m/K)                            
-              lambda_w, &   !< hydraulic diffusivity of soil (?)
-              gamma_w,  &   !< hydraulic conductivity of soil (W/m/K)
-              rho_c_total   !< volumetric heat capacity of the actual soil matrix (?) 
-
-#if defined( __nopointer )
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
-              t_soil,    & !< Soil temperature (K)
-              t_soil_av, & !< Average of t_soil
-              t_soil_p,  & !< Prog. soil temperature (K)
-              m_soil,    & !< Soil moisture (m3/m3)
-              m_soil_av, & !< Average of m_soil
-              m_soil_p     !< Prog. soil moisture (m3/m3)
-#else
-    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
-              t_soil, t_soil_p, &
-              m_soil, m_soil_p    
-
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
-              t_soil_av, t_soil_1, t_soil_2,                                   &
-              m_soil_av, m_soil_1, m_soil_2
-#endif
-
-
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
-              tt_soil_m, & !< t_soil storage array
-              tm_soil_m, & !< m_soil storage array
-              root_fr      !< root fraction (sum=1) 
 
 
@@ -571,5 +593,5 @@
 !
 !-- Public parameters, constants and initial values
-    PUBLIC land_surface, skip_time_do_lsm
+    PUBLIC aero_resist_kray, skip_time_do_lsm
 
 !
@@ -578,11 +600,6 @@
 
 !
-!-- Public 2D output variables
-    PUBLIC ghf_eb, qsws_eb, qsws_liq_eb, qsws_soil_eb,qsws_veg_eb, r_a, r_s,   &
-           shf_eb
-
-!
 !-- Public prognostic variables
-    PUBLIC m_soil, t_soil
+    PUBLIC m_soil_h, t_soil_h
 
 
@@ -871,5 +888,5 @@
     USE control_parameters,                                                    &
         ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
-               most_method, topography
+               most_method
                  
     USE radiation_model_mod,                                                   &
@@ -879,5 +896,5 @@
     IMPLICIT NONE
 
- 
+    INTEGER(iwp) ::  k        !< running index, z-dimension
 !
 !-- Dirichlet boundary conditions are required as the surface fluxes are
@@ -897,10 +914,4 @@
     ENDIF
 
-    IF ( topography /= 'flat' )  THEN
-       message_string = 'lsm cannot be used ' //                            & 
-                        'in combination with  topography /= "flat"'
-       CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
-    ENDIF
-
     IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
            most_method == 'lookup' )  THEN
@@ -1044,4 +1055,14 @@
 
     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
@@ -1059,5 +1080,5 @@
 !> Solver for the energy balance at the surface.
 !------------------------------------------------------------------------------!
- SUBROUTINE lsm_energy_balance
+ SUBROUTINE lsm_energy_balance( horizontal, l )
 
 
@@ -1065,6 +1086,15 @@
 
     INTEGER(iwp) ::  i         !< running index
+    INTEGER(iwp) ::  i_off     !< offset to determine index of surface element, seen from atmospheric grid point, for x
     INTEGER(iwp) ::  j         !< running index
-    INTEGER(iwp) ::  k, ks     !< running index
+    INTEGER(iwp) ::  j_off     !< offset to determine index of surface element, seen from atmospheric grid point, for y
+    INTEGER(iwp) ::  k         !< running index
+    INTEGER(iwp) ::  k_off     !< offset to determine index of surface element, seen from atmospheric grid point, for z
+    INTEGER(iwp) ::  k_rad     !< index to access radiation array
+    INTEGER(iwp) ::  ks        !< running index
+    INTEGER(iwp) ::  l         !< surface-facing index
+    INTEGER(iwp) ::  m         !< running index concerning wall elements
+
+    LOGICAL      ::  horizontal !< Flag indicating horizontal or vertical surfaces
 
     REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
@@ -1090,166 +1120,235 @@
                 qv1            !< specific humidity at first grid level
 
+    TYPE(surf_type_lsm), POINTER ::  surf_t_surface
+    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_soil
+    TYPE(surf_type_lsm), POINTER ::  surf_t_soil
+
+    TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable
+
+    IF ( horizontal )  THEN
+       surf              => surf_lsm_h
+
+       surf_t_surface    => t_surface_h
+       surf_t_surface_p  => t_surface_h_p
+       surf_tt_surface_m => tt_surface_h_m
+       surf_m_liq_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_soil       => m_soil_h
+       surf_t_soil       => t_soil_h
+
+       k_off     = -1
+       j_off     = 0
+       i_off     = 0
+    ELSE
+       surf              => surf_lsm_v(l)
+
+       surf_t_surface    => t_surface_v(l)
+       surf_t_surface_p  => t_surface_v_p(l)
+       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_soil       => m_soil_v(l)
+       surf_t_soil       => t_soil_v(l)
+
+       k_off = 0
+       IF ( l == 0 )  THEN 
+          j_off = -1 
+          i_off = 0
+       ELSEIF ( l == 1 )  THEN 
+          j_off = 1 
+          i_off = 0
+       ELSEIF ( l == 2 )  THEN 
+          j_off = 0 
+          i_off = -1
+       ELSEIF ( l == 3 )  THEN 
+          j_off = 0 
+          i_off = 1
+       ENDIF
+    ENDIF
+
 !
 !-- Calculate the exner function for the current time step
     exn = ( surface_pressure / 1000.0_wp )**0.286_wp
 
-    DO  i = nxlg, nxrg
-       DO  j = nysg, nyng
-          k = nzb_s_inner(j,i)
-
-!
-!--       Set lambda_surface according to stratification between skin layer and soil
-          IF (  .NOT.  pave_surface(j,i) )  THEN
-
-             c_surface_tmp = c_surface
-
-             IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
-                lambda_surface = lambda_surface_s(j,i)
-             ELSE
-                lambda_surface = lambda_surface_u(j,i)
-             ENDIF
+    DO  m = 1, surf%ns
+
+       i   = surf%i(m)            
+       j   = surf%j(m)
+       k   = surf%k(m)
+!
+!--    Determine height index for radiation. Note, in clear-sky case radiation
+!--    arrays have rank 0 in first dimensions, so index must be zero. In case
+!--    of RRTMG radiation arrays have non-zero rank in first dimension, so that 
+!--    radiation can be obtained at respective height level
+       k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )
+!
+!--    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
-
-             c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
-             lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
-
-          ENDIF
-
-!
-!--       First step: calculate aerodyamic resistance. As pt, us, ts
-!--       are not available for the prognostic time step, data from the last
-!--       time step is used here. Note that this formulation is the
-!--       equivalent to the ECMWF formulation using drag coefficients
-          IF ( cloud_physics )  THEN
-             pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
-             qv1 = q(k+1,j,i) - ql(k+1,j,i)
-          ELSE
-             pt1 = pt(k+1,j,i)
-             qv1 = q(k+1,j,i)
-          ENDIF
-
-          r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
-
-!
-!--       Make sure that the resistance does not drop to zero for neutral 
-!--       stratification
-          IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
-
-!
-!--       Second step: calculate canopy resistance r_canopy
-!--       f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
+             lambda_surface = surf%lambda_surface_u(m)
+          ENDIF
+       ELSE
+
+          c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
+          lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
+
+       ENDIF
+
+!
+!--    First step: calculate aerodyamic resistance. As pt, us, ts
+!--    are not available for the prognostic time step, data from the last
+!--    time step is used here. Note that this formulation is the
+!--    equivalent to the ECMWF formulation using drag coefficients
+       IF ( cloud_physics )  THEN
+          pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+          qv1 = q(k,j,i) - ql(k,j,i)
+       ELSE
+          pt1 = pt(k,j,i)
+          qv1 = q(k,j,i)
+       ENDIF
+!
+!--    Calculate aerodynamical resistance. For horizontal and vertical 
+!--    surfaces MOST is applied. Moreover, for vertical surfaces, resistance
+!--    can be obtain via parameterization of Mason (2000) / 
+!--    Krayenhoff and Voogt (2006).
+!--    To do: detailed investigation which approach is better!
+       IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
+          surf%r_a(m) = ( pt1 - surf_lsm_h%pt_surface(m) ) /                   &
+                        ( surf%ts(m) * surf%us(m) + 1.0E-20_wp )
+       ELSE
+          surf%r_a(m) = 1.0_wp / ( 11.8_wp + 4.2_wp *                          &
+                        SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +      &
+                              ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +      &
+                              ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 )      &
+                                 )
+       ENDIF
+!
+!--    Make sure that the resistance does not drop to zero for neutral 
+!--    stratification
+       IF ( ABS( surf%r_a(m) ) < 1.0_wp )  surf%r_a(m) = 1.0_wp
+!
+!--    Second step: calculate canopy resistance r_canopy
+!--    f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
  
-!--       f1: correction for incoming shortwave radiation (stomata close at 
-!--       night)
-          f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /        &
-                           (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)             &
-                            + 1.0_wp)) )
-
-
-
-!
-!--       f2: correction for soil moisture availability to plants (the 
-!--       integrated soil moisture must thus be considered here)
-!--       f2 = 0 for very dry soils
-          m_total = 0.0_wp
-          DO  ks = nzb_soil, nzt_soil
-              m_total = m_total + root_fr(ks,j,i)                              &
-                        * MAX(m_soil(ks,j,i),m_wilt(j,i))
-          ENDDO 
-
-          IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
-             f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
-          ELSEIF ( m_total >= m_fc(j,i) )  THEN
-             f2 = 1.0_wp
-          ELSE
-             f2 = 1.0E-20_wp
-          ENDIF
-
-!
-!--       Calculate water vapour pressure at saturation 
-          e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)        &
-                        - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
-
-!
-!--       f3: correction for vapour pressure deficit
-          IF ( g_d(j,i) /= 0.0_wp )  THEN
-!
-!--          Calculate vapour pressure
-             e  = qv1 * surface_pressure / 0.622_wp
-             f3 = EXP ( -g_d(j,i) * (e_s - e) )
-          ELSE
-             f3 = 1.0_wp
-          ENDIF
-
-!
-!--       Calculate canopy resistance. In case that c_veg is 0 (bare soils),
-!--       this calculation is obsolete, as r_canopy is not used below.
-!--       To do: check for very dry soil -> r_canopy goes to infinity
-          r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * 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 ( soil_type_2d(j,i) /= 7 )  THEN
-             m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *        &
-                     m_res(j,i)
-          ELSE
-             m_min = m_wilt(j,i)
-          ENDIF
-
-          f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
-          f2 = MAX(f2,1.0E-20_wp)
-          f2 = MIN(f2,1.0_wp)
-
-          r_soil(j,i) = r_soil_min(j,i) / f2
-
-!
-!--       Calculate the maximum possible liquid water amount on plants and
-!--       bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
-!--       assumed, while paved surfaces might hold up 1 mm of water. The 
-!--       liquid water fraction for paved surfaces is calculated after 
-!--       Noilhan & Planton (1989), while the ECMWF formulation is used for
-!--       vegetated surfaces and bare soils.
-          IF ( pave_surface(j,i) )  THEN
-             m_liq_eb_max = m_max_depth * 5.0_wp
-             c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
-          ELSE
-             m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)              &
-                            + (1.0_wp - c_veg(j,i)) )
-             c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
-          ENDIF
-
-!
-!--       Calculate saturation specific humidity
-          q_s = 0.622_wp * e_s / surface_pressure
-
-!
-!--       In case of dewfall, set evapotranspiration to zero
-!--       All super-saturated water is then removed from the air
-          IF ( humidity  .AND.  q_s <= qv1 )  THEN
-             r_canopy(j,i) = 0.0_wp
-             r_soil(j,i)   = 0.0_wp
-          ENDIF
-
-!
-!--       Calculate coefficients for the total evapotranspiration 
-!--       In case of water surface, set vegetation and soil fluxes to zero.
-!--       For pavements, only evaporation of liquid water is possible.
-          IF ( water_surface(j,i) )  THEN
-             f_qsws_veg  = 0.0_wp
-             f_qsws_soil = 0.0_wp
-             f_qsws_liq  = rho_lv / r_a(j,i)
-          ELSEIF ( pave_surface (j,i) )  THEN
-             f_qsws_veg  = 0.0_wp
-             f_qsws_soil = 0.0_wp
-             f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
-          ELSE
-             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
-                           (r_a(j,i) + r_canopy(j,i))
-             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
-                                                             r_soil(j,i))
-             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
-          ENDIF
+!--    f1: correction for incoming shortwave radiation (stomata close at 
+!--    night)
+       f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off) + 0.05_wp ) /&
+                        (0.81_wp * (0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off)&
+                         + 1.0_wp)) )
+!
+!--    f2: correction for soil moisture availability to plants (the 
+!--    integrated soil moisture must thus be considered here)
+!--    f2 = 0 for very dry soils
+       m_total = 0.0_wp
+       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) )
+       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
+          f2 = 1.0_wp
+       ELSE
+          f2 = 1.0E-20_wp
+       ENDIF
+!
+!--    Calculate water vapour pressure at saturation 
+       e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( surf_t_surface%var_1d(m) &
+                     - 273.16_wp ) / ( surf_t_surface%var_1d(m) - 35.86_wp ) )
+!
+!--    f3: correction for vapour pressure deficit
+       IF ( surf%g_d(m) /= 0.0_wp )  THEN
+!
+!--       Calculate vapour pressure
+          e  = qv1 * surface_pressure / 0.622_wp
+          f3 = EXP ( - surf%g_d(m) * (e_s - e) )
+       ELSE
+          f3 = 1.0_wp
+       ENDIF
+!
+!--    Calculate canopy resistance. In case that c_veg is 0 (bare soils),
+!--    this calculation is obsolete, as r_canopy is not used below.
+!--    To do: check for very dry soil -> r_canopy goes to infinity
+       surf%r_canopy(m) = surf%r_canopy_min(m) /                               &
+                              ( 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 )
+       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
+!--    bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
+!--    assumed, while paved surfaces might hold up 1 mm of water. The 
+!--    liquid water fraction for paved surfaces is calculated after 
+!--    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 )
+       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 )
+       ENDIF
+!
+!--    Calculate saturation specific humidity
+       q_s = 0.622_wp * e_s / surface_pressure
+!
+!--    In case of dewfall, set evapotranspiration to zero
+!--    All super-saturated water is then removed from the air
+       IF ( humidity  .AND.  q_s <= qv1 )  THEN
+          surf%r_canopy(m) = 0.0_wp
+          surf%r_soil(m)   = 0.0_wp
+       ENDIF
+
+!
+!--    Calculate coefficients for the total evapotranspiration 
+!--    In case of water surface, set vegetation and soil fluxes to zero.
+!--    For pavements, only evaporation of liquid water is possible.
+       IF ( surf%water_surface(m) )  THEN
+          f_qsws_veg  = 0.0_wp
+          f_qsws_soil = 0.0_wp
+          f_qsws_liq  = rho_lv / surf%r_a(m)
+       ELSEIF ( surf%pave_surface (m) )  THEN
+          f_qsws_veg  = 0.0_wp
+          f_qsws_soil = 0.0_wp
+          f_qsws_liq  = rho_lv * surf%c_liq(m) / surf%r_a(m)
+       ELSE
+          f_qsws_veg  = rho_lv * surf%c_veg(m) *                               &
+                            ( 1.0_wp        - surf%c_liq(m)    ) /             &
+                            ( surf%r_a(m) + surf%r_canopy(m) )
+          f_qsws_soil = rho_lv * (1.0_wp    - surf%c_veg(m)    ) /             &
+                            ( surf%r_a(m) + surf%r_soil(m)   )
+          f_qsws_liq  = rho_lv * surf%c_veg(m) * surf%c_liq(m)   /             &
+                              surf%r_a(m)
+       ENDIF
 !
 !--       If soil moisture is below wilting point, plants do no longer
@@ -1259,233 +1358,228 @@
 !           ENDIF
 
-          f_shf  = rho_cp / r_a(j,i)
-          f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
-
-!
-!--       Calculate derivative of q_s for Taylor series expansion
-          e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -           &
-                           17.269_wp*(t_surface(j,i) - 273.16_wp)              &
-                           / (t_surface(j,i) - 35.86_wp)**2 )
-
-          dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
-
-!
-!--       Add LW up so that it can be removed in prognostic equation
-          rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
-
-!
-!--       Calculate new skin temperature
-          IF ( humidity )  THEN
-
-!
-!--          Numerator of the prognostic equation
-             coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
-                      * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
-                      + f_shf * pt1 + f_qsws * ( qv1 - q_s                     &
-                      + dq_s_dt * t_surface(j,i) ) + lambda_surface            &
-                      * t_soil(nzb_soil,j,i)
-
-!
-!--          Denominator of the prognostic equation
-             coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt              &
-                      + lambda_surface + f_shf / exn
-          ELSE
-
-!
-!--          Numerator of the prognostic equation
-             coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
-                      * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
-                      + f_shf * pt1  + lambda_surface                          &
-                      * t_soil(nzb_soil,j,i)
-
-!
-!--          Denominator of the prognostic equation
-             coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
-
-          ENDIF
-
-          tend = 0.0_wp
-
-!
-!--       Implicit solution when the surface layer has no heat capacity,
-!--       otherwise use RK3 scheme.
-          t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *       &
-                             t_surface(j,i) ) / ( c_surface_tmp + coef_2       &
-                                * dt_3d * tsc(2) ) 
-
-!
-!--       Add RK3 term
-          IF ( c_surface_tmp /= 0.0_wp )  THEN
-
-             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
-                                * tt_surface_m(j,i)
-
-!
-!--          Calculate true tendency
-             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
-                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
-!
-!--          Calculate t_surface tendencies for the next Runge-Kutta step
-             IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
-                   tt_surface_m(j,i) = tend
-                ELSEIF ( intermediate_timestep_count <                         &
-                         intermediate_timestep_count_max )  THEN
-                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
-                                       * tt_surface_m(j,i)
-                ENDIF
+       f_shf  = rho_cp / surf%r_a(m)
+       f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
+!
+!--    Calculate derivative of q_s for Taylor series expansion
+       e_s_dt = e_s * ( 17.269_wp / ( surf_t_surface%var_1d(m) - 35.86_wp) -   &
+                        17.269_wp*( surf_t_surface%var_1d(m) - 273.16_wp)      &
+                       / ( surf_t_surface%var_1d(m) - 35.86_wp)**2 )
+
+       dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
+!
+!--    Add LW up so that it can be removed in prognostic equation
+       surf%rad_net_l(m) = rad_net(j,i) + rad_lw_out(nzb,j,i)
+!
+!--    Calculate new skin temperature
+       IF ( humidity )  THEN
+!
+!--       Numerator of the prognostic equation
+          coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
+                   * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
+                   + f_shf * pt1 + f_qsws * ( qv1 - q_s                        &
+                   + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface     &
+                   * surf_t_soil%var_2d(nzb_soil,m)
+!
+!--       Denominator of the prognostic equation
+          coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt                 &
+                   + lambda_surface + f_shf / exn
+       ELSE
+!
+!--       Numerator of the prognostic equation
+          coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
+                   * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
+                   + f_shf * pt1  + lambda_surface                             &
+                   * surf_t_soil%var_2d(nzb_soil,m)
+!
+!--       Denominator of the prognostic equation
+          coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
+
+       ENDIF
+
+       tend = 0.0_wp
+
+!
+!--    Implicit solution when the surface layer has no heat capacity,
+!--    otherwise use RK3 scheme.
+       surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
+                          surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
+                                             * dt_3d * tsc(2) ) 
+
+!
+!--    Add RK3 term
+       IF ( c_surface_tmp /= 0.0_wp )  THEN
+
+          surf_t_surface_p%var_1d(m) = surf_t_surface_p%var_1d(m) + dt_3d *    &
+                                       tsc(3) * surf_tt_surface_m%var_1d(m)
+
+!
+!--       Calculate true tendency
+          tend = ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) -     &
+                   dt_3d * tsc(3) * surf_tt_surface_m%var_1d(m)) / (dt_3d  * tsc(2))
+!
+!--       Calculate t_surface tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                surf_tt_surface_m%var_1d(m) = tend
+             ELSEIF ( intermediate_timestep_count <                            &
+                      intermediate_timestep_count_max )  THEN
+                surf_tt_surface_m%var_1d(m) = -9.5625_wp * tend +              &
+                                               5.3125_wp * surf_tt_surface_m%var_1d(m)
              ENDIF
           ENDIF
-
-!
-!--       In case of fast changes in the skin temperature, it is possible to
-!--       update the radiative fluxes independently from the prescribed
-!--       radiation call frequency. This effectively prevents oscillations,
-!--       especially when setting skip_time_do_radiation /= 0. The threshold
-!--       value of 0.2 used here is just a first guess. This method should be
-!--       revised in the future as tests have shown that the threshold is 
-!--       often reached, when no oscillations would occur (causes immense
-!--       computing time for the radiation code).
-          IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.        &
-               unscheduled_radiation_calls )  THEN
-             force_radiation_call_l = .TRUE.
-          ENDIF
-
-          pt(k,j,i) = t_surface_p(j,i) / exn
-
-!
-!--       Calculate fluxes
-          rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)         &
-                             * t_surface(j,i) - rad_lw_out(nzb,j,i)            &
-                             - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
-
-          rad_net(j,i) = rad_net_l(j,i)
-          rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) &
-                                * ( t_surface_p(j,i) - t_surface(j,i) )
-
-          ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)                  &
-                           - t_soil(nzb_soil,j,i))
-
-          shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
-
-          shf(j,i) = shf_eb(j,i) / rho_cp
-
-          IF ( humidity )  THEN
-             qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt               &
-                             * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
-
-             qsws(j,i) = qsws_eb(j,i) / rho_lv
-
-             qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                   &
-                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
-                                 * t_surface_p(j,i) )
-
-             qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                   &
-                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
-                                 * t_surface_p(j,i) )
-
-             qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                   &
-                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
-                                 * t_surface_p(j,i) )
-          ENDIF
-
-!
-!--       Calculate the true surface resistance
-          IF ( qsws_eb(j,i) == 0.0_wp )  THEN
-             r_s(j,i) = 1.0E10_wp
-          ELSE
-             r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
-                        * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )        &
-                        / qsws_eb(j,i) - r_a(j,i)
-          ENDIF
-
-!
-!--       Calculate change in liquid water reservoir due to dew fall or 
-!--       evaporation of liquid water
-          IF ( humidity )  THEN
-!
-!--          If precipitation is activated, add rain water to qsws_liq_eb
-!--          and qsws_soil_eb according the the vegetation coverage.
-!--          precipitation_rate is given in mm.
-             IF ( precipitation )  THEN
-
-!
-!--             Add precipitation to liquid water reservoir, if possible.
-!--             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 ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
-                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                         &
-                                    + c_veg(j,i) * prr(k,j,i) * hyrho(k)       &
-                                    * 0.001_wp * rho_l * l_v
-                ELSE
-                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
-                                     + c_veg(j,i) * prr(k,j,i) * hyrho(k)      &
-                                     * 0.001_wp * rho_l * l_v
-                ENDIF
-
-!--             Add precipitation to bare soil according to the bare soil
-!--             coverage.
-                qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp                &
-                                    - c_veg(j,i)) * prr(k,j,i) * hyrho(k)      &
-                                    * 0.001_wp * rho_l * l_v
+       ENDIF
+
+!
+!--    In case of fast changes in the skin temperature, it is possible to
+!--    update the radiative fluxes independently from the prescribed
+!--    radiation call frequency. This effectively prevents oscillations,
+!--    especially when setting skip_time_do_radiation /= 0. The threshold
+!--    value of 0.2 used here is just a first guess. This method should be
+!--    revised in the future as tests have shown that the threshold is 
+!--    often reached, when no oscillations would occur (causes immense
+!--    computing time for the radiation code).
+       IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) ) > 0.2_wp  .AND. &
+            unscheduled_radiation_calls )  THEN
+          force_radiation_call_l = .TRUE.
+       ENDIF
+
+       pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
+       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exn
+
+!
+!--    Calculate fluxes
+       surf%rad_net_l(m) = surf%rad_net_l(m) +                                 &
+                            rad_lw_out_change_0(j,i)                           &
+                          * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)     &
+                          - rad_lw_out_change_0(j,i) * surf_t_surface_p%var_1d(m)
+
+       rad_net(j,i) = surf%rad_net_l(m)
+       rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) *  &
+                     ( surf_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_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
+
+       IF ( humidity )  THEN
+          surf%qsws_eb(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                   &
+                              + 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                   &
+                              + 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                   &
+                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
+                              * surf_t_surface_p%var_1d(m) )
+       ENDIF
+
+!
+!--    Calculate the true surface resistance
+       IF ( surf%qsws_eb(m) == 0.0_wp )  THEN
+          surf%r_s(m) = 1.0E10_wp
+       ELSE
+          surf%r_s(m) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
+                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
+                            surf_t_surface_p%var_1d(m) ) /                     &
+                            surf%qsws_eb(m) - surf%r_a(m)
+       ENDIF
+
+!
+!--    Calculate change in liquid water reservoir due to dew fall or 
+!--    evaporation of liquid water
+       IF ( humidity )  THEN
+!
+!--       If precipitation is activated, add rain water to qsws_liq_eb
+!--       and qsws_soil_eb according the the vegetation coverage.
+!--       precipitation_rate is given in mm.
+          IF ( precipitation )  THEN
+
+!
+!--          Add precipitation to liquid water reservoir, if possible.
+!--          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)                      &
+                                 + 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%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
+                                 * hyrho(k+k_off)                              &
+                                 * 0.001_wp * rho_l * l_v
              ENDIF
 
-!
-!--          If the air is saturated, check the reservoir water level
-             IF ( qsws_eb(j,i) < 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),       
-!--             so that tend is zero and no further check is needed
-                IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
-                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
-                                        + qsws_liq_eb(j,i)
-
-                   qsws_liq_eb(j,i)  = 0.0_wp
-                ENDIF
-
-!
-!--             In case qsws_veg_eb becomes negative (unphysical behavior), 
-!--             let the water enter the liquid water reservoir as dew on the
-!--             plant
-                IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
-                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
-                   qsws_veg_eb(j,i) = 0.0_wp
-                ENDIF
-             ENDIF                    
-  
-             tend = - qsws_liq_eb(j,i) * drho_l_lv
-             m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend         &
-                                                + tsc(3) * tm_liq_eb_m(j,i) )
-
-!
-!--          Check if reservoir is overfull -> reduce to maximum
-!--          (conservation of water is violated here)
-             m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
-
-!
-!--          Check if reservoir is empty (avoid values < 0.0)
-!--          (conservation of water is violated here)
-             m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
-
-
-!
-!--          Calculate m_liq_eb tendencies for the next Runge-Kutta step
-             IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
-                   tm_liq_eb_m(j,i) = tend
-                ELSEIF ( intermediate_timestep_count <                         &
-                         intermediate_timestep_count_max )  THEN
-                   tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
-                                    * tm_liq_eb_m(j,i)
-                ENDIF
+!--          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%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
+                               * hyrho(k+k_off)                                &
+                               * 0.001_wp * rho_l * l_v
+          ENDIF
+
+!
+!--       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),       
+!--          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
              ENDIF
 
-          ENDIF
-
-       ENDDO
+!
+!--          In case qsws_veg_eb 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
+             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 *        &
+                                        ( tsc(2) * tend +                      &
+                                          tsc(3) * surf_tm_liq_eb_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 )
+
+!
+!--       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
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                surf_tm_liq_eb_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)
+             ENDIF
+          ENDIF
+
+       ENDIF
+
     ENDDO
 
@@ -1497,5 +1591,5 @@
 #if defined( __parallel )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
+       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,       &
                            1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
 #else
@@ -1513,5 +1607,59 @@
 !
 !-- Calculate new roughness lengths (for water surfaces only)
-    CALL calc_z0_water_surface
+    IF ( horizontal )  CALL calc_z0_water_surface
+
+    CONTAINS
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculation of specific humidity of the skin layer (surface). It is assumend
+!> that the skin is always saturated.
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_q_surface
+
+       IMPLICIT NONE
+
+       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
+
+       DO  m = 1, surf%ns
+
+          i   = surf%i(m)            
+          j   = surf%j(m)
+          k   = surf%k(m)
+
+!
+!--       Calculate water vapour pressure at saturation
+          e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp *                         &
+                                 ( surf_t_surface_p%var_1d(m) - 273.16_wp ) /  &
+                                 ( surf_t_surface_p%var_1d(m) - 35.86_wp  )    &
+                                         )
+
+!
+!--       Calculate specific humidity at saturation
+          q_s = 0.622_wp * e_s / surface_pressure
+
+          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )
+
+!
+!--       Calculate specific humidity at surface
+          IF ( cloud_physics )  THEN
+             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
+                                        ( 1.0_wp - resistance ) *              &
+                                        ( q(k,j,i) - ql(k,j,i) )
+          ELSE
+             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
+                                        ( 1.0_wp - resistance ) *              &
+                                          q(k,j,i)
+          ENDIF
+
+!
+!--       Update virtual potential temperature
+          vpt(k+k_off,j+j_off,i+i_off) = pt(k+k_off,j+j_off,i+i_off) *         &
+                     ( 1.0_wp + 0.61_wp * q(k+k_off,j+j_off,i+i_off) )
+
+       ENDDO
+
+    END SUBROUTINE calc_q_surface
+
 
 
@@ -1607,16 +1755,17 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i !< running index
-       INTEGER(iwp) ::  j !< running index
-       INTEGER(iwp) ::  k !< running index
+       INTEGER(iwp) ::  i     !< running index
+       INTEGER(iwp) ::  i_off !< index offset of surface element, seen from atmospheric grid point 
+       INTEGER(iwp) ::  j     !< running index
+       INTEGER(iwp) ::  j_off !< index offset of surface element, seen from atmospheric grid point 
+       INTEGER(iwp) ::  k     !< running index
+       INTEGER(iwp) ::  l     !< running index surface facing
+       INTEGER(iwp) ::  m     !< running index
 
        REAL(wp) :: pt1   !< potential temperature at first grid level
-
 
 !
 !--    Calculate Exner function
        exn = ( surface_pressure / 1000.0_wp )**0.286_wp
-
-
 !
 !--    If no cloud physics is used, rho_surface has not been calculated before
@@ -1634,43 +1783,89 @@
 !
 !--    Set inital values for prognostic quantities
-       tt_surface_m = 0.0_wp
-       tt_soil_m    = 0.0_wp
-       tm_soil_m    = 0.0_wp
-       tm_liq_eb_m  = 0.0_wp
-       c_liq        = 0.0_wp
-
-       ghf_eb = 0.0_wp
-       shf_eb = rho_cp * shf
+!--    Horizontal surfaces
+       tt_surface_h_m%var_1d = 0.0_wp
+       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
+       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
-          qsws_eb = rho_lv * qsws
+          surf_lsm_h%qsws_eb = rho_lv * surf_lsm_h%qsws
        ELSE
-          qsws_eb = 0.0_wp
-       ENDIF
-
-       qsws_liq_eb  = 0.0_wp
-       qsws_soil_eb = 0.0_wp
-       qsws_veg_eb  = 0.0_wp
-
-       r_a        = 50.0_wp
-       r_s        = 50.0_wp
-       r_canopy   = 0.0_wp
-       r_soil     = 0.0_wp
+          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%r_a        = 50.0_wp
+       surf_lsm_h%r_s        = 50.0_wp
+       surf_lsm_h%r_canopy   = 0.0_wp
+       surf_lsm_h%r_soil     = 0.0_wp
+!
+!--    Do the same for vertical surfaces
+       DO  l = 0, 3
+          tt_surface_v_m(l)%var_1d = 0.0_wp
+          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
+          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)%r_a        = 50.0_wp
+          surf_lsm_v(l)%r_s        = 50.0_wp
+          surf_lsm_v(l)%r_canopy   = 0.0_wp
+          surf_lsm_v(l)%r_soil     = 0.0_wp
+       ENDDO
 
 !
 !--    Allocate 3D soil model arrays
-       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-
-       lambda_h = 0.0_wp
+!--    First, for horizontal surfaces
+       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
 !
 !--    If required, allocate humidity-related variables for the soil model
        IF ( humidity )  THEN
-          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
-
-          lambda_w = 0.0_wp 
-       ENDIF
+          ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
+          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )  
+
+          surf_lsm_h%lambda_w = 0.0_wp 
+       ENDIF
+!
+!--    For vertical surfaces
+       DO  l = 0, 3
+          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
+!
+!--       If required, allocate humidity-related variables for the soil model
+          IF ( humidity )  THEN
+             ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
+             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )  
+
+             surf_lsm_v(l)%lambda_w = 0.0_wp 
+          ENDIF      
+       ENDDO
 
 !
@@ -1735,21 +1930,38 @@
 !
 !--    Map values to the respective 2D arrays
-       alpha_vg      = alpha_vangenuchten
-       l_vg          = l_vangenuchten
-       n_vg          = n_vangenuchten 
-       gamma_w_sat   = hydraulic_conductivity
-       m_sat         = saturation_moisture
-       m_fc          = field_capacity
-       m_wilt        = wilting_point
-       m_res         = residual_moisture
-       r_soil_min    = min_soil_resistance
+!--    Horizontal surfaces
+       surf_lsm_h%alpha_vg      = alpha_vangenuchten
+       surf_lsm_h%l_vg          = l_vangenuchten
+       surf_lsm_h%n_vg          = n_vangenuchten 
+       surf_lsm_h%gamma_w_sat   = hydraulic_conductivity
+       surf_lsm_h%m_sat         = saturation_moisture
+       surf_lsm_h%m_fc          = field_capacity
+       surf_lsm_h%m_wilt        = wilting_point
+       surf_lsm_h%m_res         = residual_moisture
+       surf_lsm_h%r_soil_min    = min_soil_resistance
+!
+!--    Vertical surfaces
+       DO  l = 0, 3
+          surf_lsm_v(l)%alpha_vg      = alpha_vangenuchten
+          surf_lsm_v(l)%l_vg          = l_vangenuchten
+          surf_lsm_v(l)%n_vg          = n_vangenuchten 
+          surf_lsm_v(l)%gamma_w_sat   = hydraulic_conductivity
+          surf_lsm_v(l)%m_sat         = saturation_moisture
+          surf_lsm_v(l)%m_fc          = field_capacity
+          surf_lsm_v(l)%m_wilt        = wilting_point
+          surf_lsm_v(l)%m_res         = residual_moisture
+          surf_lsm_v(l)%r_soil_min    = min_soil_resistance
+       ENDDO
+
+
 
 !
 !--    Initial run actions
        IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
-
-          t_soil    = 0.0_wp
-          m_liq_eb  = 0.0_wp
-          m_soil    = 0.0_wp
+!
+!--       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
 
 !
@@ -1758,34 +1970,82 @@
 !--       wilting point) -> problems with devision by zero)
           DO  k = nzb_soil, nzt_soil
-             t_soil(k,:,:)    = soil_temperature(k)
-             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
+             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(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
+          t_soil_h%var_2d(nzt_soil+1,:) = soil_temperature(nzt_soil+1)
 
 !
 !--       Calculate surface temperature
-          t_surface   = pt_surface * exn
+          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
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                k = nzb_s_inner(j,i)
+!--       for the first time step. Only for interior core domain, not for ghost points
+          DO  m = 1, surf_lsm_h%ns
+             i   = surf_lsm_h%i(m)            
+             j   = surf_lsm_h%j(m)
+             k   = surf_lsm_h%k(m)
+
+             IF ( cloud_physics )  THEN
+                pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+             ELSE
+                pt1 = pt(k,j,i)
+             ENDIF
+
+!
+!--          Assure that r_a cannot be zero at model start
+             IF ( pt1 == surf_lsm_h%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
+
+             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
+             DO  m = 1, surf_lsm_v(l)%ns
+                i   = surf_lsm_v(l)%i(m)            
+                j   = surf_lsm_v(l)%j(m)
+                k   = surf_lsm_v(l)%k(m)
 
                 IF ( cloud_physics )  THEN
-                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+                   pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
                 ELSE
-                   pt1 = pt(k+1,j,i)
+                   pt1 = pt(k,j,i)
                 ENDIF
 
 !
 !--             Assure that r_a cannot be zero at model start
-                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
-
-                us(j,i)  = 0.1_wp
-                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
-                shf(j,i) = - us(j,i) * ts(j,i)
+                IF ( pt1 == surf_lsm_h%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)
              ENDDO
           ENDDO
@@ -1794,16 +2054,61 @@
 !--    Actions for restart runs
        ELSE
-
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                k = nzb_s_inner(j,i)                
-                t_surface(j,i) = pt(k,j,i) * exn
-             ENDDO
+!
+!--       Horizontal surfaces
+          DO  m = 1, surf_lsm_h%ns
+             i   = surf_lsm_h%i(m)            
+             j   = surf_lsm_h%j(m)
+             k   = surf_lsm_h%k(m)          
+             t_surface_h%var_1d(m) = pt(k-1,j,i) * exn
           ENDDO
-
-       ENDIF
-
-       DO  k = nzb_soil, nzt_soil
-          root_fr(k,:,:) = root_fraction(k)
+!
+!--       Vertical surfaces
+          DO  l = 0, 3
+!
+!--          Set index offset of surface element, seen from atmospheric grid point 
+             IF ( l == 0 )  THEN 
+                j_off = -1 
+                i_off = 0
+             ELSEIF ( l == 1 )  THEN 
+                j_off = 1 
+                i_off = 0
+             ELSEIF ( l == 2 )  THEN 
+                j_off = 0 
+                i_off = -1
+             ELSEIF ( l == 3 )  THEN 
+                j_off = 0 
+                i_off = 1
+             ENDIF
+             DO  m = 1, surf_lsm_v(l)%ns
+                i   = surf_lsm_v(l)%i(m)            
+                j   = surf_lsm_v(l)%j(m)
+                k   = surf_lsm_v(l)%k(m)          
+                t_surface_v(l)%var_1d(m) = pt(k,j+j_off,i+i_off) * exn
+             ENDDO
+          ENDDO
+
+       ENDIF
+!
+!--    Initialize root fraction
+!--    Horizontal surfaces
+       DO  m = 1, surf_lsm_h%ns
+          i   = surf_lsm_h%i(m)            
+          j   = surf_lsm_h%j(m)
+
+          DO  k = nzb_soil, nzt_soil
+             surf_lsm_h%root_fr(k,m) = root_fraction(k)
+          ENDDO
+       ENDDO
+!
+!--    Vertical surfaces
+       DO  l = 0, 3
+          DO  m = 1, surf_lsm_v(l)%ns
+             i   = surf_lsm_v(l)%i(m)            
+             j   = surf_lsm_v(l)%j(m)
+
+             DO  k = nzb_soil, nzt_soil
+                surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
+             ENDDO
+          ENDDO
        ENDDO
 
@@ -1843,7 +2148,23 @@
 
           IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
-             DO  k = nzb_soil, nzt_soil
-                root_fr(k,:,:) = root_distribution(k,veg_type)
-                root_fraction(k) = root_distribution(k,veg_type)
+             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)
+                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
@@ -1877,19 +2198,39 @@
 !--    Map vegetation and soil types to 2D array to allow for heterogeneous 
 !--    surfaces via user interface see below
-       veg_type_2d = veg_type
-       soil_type_2d = soil_type
+!--    First for horizontal surfaces
+       surf_lsm_h%veg_type_2d  = veg_type
+       surf_lsm_h%soil_type_2d = soil_type
 
 !
 !--    Map vegetation parameters to the respective 2D arrays
-       r_canopy_min         = min_canopy_resistance
-       lai                  = leaf_area_index
-       c_veg                = vegetation_coverage
-       g_d                  = canopy_resistance_coefficient
-       lambda_surface_s     = lambda_surface_stable
-       lambda_surface_u     = lambda_surface_unstable
-       f_sw_in              = f_shortwave_incoming
-       z0                   = z0_eb
-       z0h                  = z0h_eb
-       z0q                  = z0q_eb
+       surf_lsm_h%r_canopy_min      = min_canopy_resistance
+       surf_lsm_h%lai               = leaf_area_index
+       surf_lsm_h%c_veg             = vegetation_coverage
+       surf_lsm_h%g_d               = canopy_resistance_coefficient
+       surf_lsm_h%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
+
+!
+!--       Map vegetation parameters to the respective 2D arrays
+          surf_lsm_v(l)%r_canopy_min      = min_canopy_resistance
+          surf_lsm_v(l)%lai               = leaf_area_index
+          surf_lsm_v(l)%c_veg             = vegetation_coverage
+          surf_lsm_v(l)%g_d               = canopy_resistance_coefficient
+          surf_lsm_v(l)%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
 
 !
@@ -1900,25 +2241,34 @@
 !--    Set flag parameter if vegetation type was set to a water surface. Also
 !--    set temperature to a constant value in all "soil" layers.
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             IF ( veg_type_2d(j,i) == 14  .OR.  veg_type_2d(j,i) == 15 )  THEN
-                water_surface(j,i) = .TRUE.
-                t_soil(:,j,i) = t_surface(j,i)
-             ELSEIF ( veg_type_2d(j,i) == 20 )  THEN
-                pave_surface(j,i) = .TRUE.
-                m_soil(:,j,i) = 0.0_wp
-             ENDIF
-
-          ENDDO
+!--    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
 
 !
-!--    Calculate new roughness lengths (for water surfaces only)
+!--    Calculate new roughness lengths (for water surfaces only, i.e. only 
+!-     horizontal surfaces)
        CALL calc_z0_water_surface
 
-       t_soil_p    = t_soil
-       m_soil_p    = m_soil
-       m_liq_eb_p  = m_liq_eb
-       t_surface_p = t_surface
+       t_soil_h_p    = t_soil_h
+       m_soil_h_p    = m_soil_h
+       m_liq_eb_h_p  = m_liq_eb_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
+       t_surface_v_p = t_surface_v
 
 
@@ -1926,10 +2276,8 @@
 !--    Store initial profiles of t_soil and m_soil (assuming they are 
 !--    horizontally homogeneous on this PE)
-       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
-                                                nysg,nxlg), 2,                 &
-                                                statistic_regions+1 )
-       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
-                                                nysg,nxlg), 2,                 &
-                                                statistic_regions+1 )
+       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
+                                                2, statistic_regions+1 )
+       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1),  & 
+                                                2, statistic_regions+1 )
 
     END SUBROUTINE lsm_init
@@ -1946,79 +2294,175 @@
        IMPLICIT NONE
 
-!
-!--    Allocate surface and soil temperature / humidity
+       INTEGER(iwp) ::  l !< index indicating facing of surface array 
+
+!
+!--    Allocate surface and soil temperature / humidity. Please note, 
+!--    these arrays are allocated according to surface-data structure,
+!--    even if they do not belong to the data type due to the 
+!--    pointer arithmetric (TARGET attribute is not allowed in a data-type). 
 #if defined( __nopointer )
-       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
+!
+!--    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 ( t_surface_h%var_1d(1:surf_lsm_h%ns)                    )
+       ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns)                  )
+       ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
+       ALLOCATE ( m_soil_h_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
+       ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns)   )
+       ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
+!
+!--    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 ( 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)                  )
+          ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
+          ALLOCATE ( m_soil_v(l)_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
+          ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns)   )
+          ALLOCATE ( t_soil_v(l)_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
+       ENDDO
 #else
-       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
+!
+!--    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 ( t_surface_h_1%var_1d(1:surf_lsm_h%ns)                  )
+       ALLOCATE ( t_surface_h_2%var_1d(1:surf_lsm_h%ns)                  )
+       ALLOCATE ( m_soil_h_1%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
+       ALLOCATE ( m_soil_h_2%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
+       ALLOCATE ( t_soil_h_1%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
+       ALLOCATE ( t_soil_h_2%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
+!
+!--    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 ( 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)                  )
+          ALLOCATE ( m_soil_v_1(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
+          ALLOCATE ( m_soil_v_2(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
+          ALLOCATE ( t_soil_v_1(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
+          ALLOCATE ( t_soil_v_2(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
+       ENDDO
 #endif
 
 !
 !--    Allocate intermediate timestep arrays
-       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+!--    Horizontal surfaces
+       ALLOCATE ( tm_liq_eb_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)  )
+       ALLOCATE ( tt_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  ) 
+!
+!--    Horizontal surfaces
+       DO  l = 0, 3
+          ALLOCATE ( tm_liq_eb_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)  )
+          ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
+       ENDDO 
+!
+!--    Allocate skin-surface temperature
+       ALLOCATE ( surf_lsm_h%pt_surface(1:surf_lsm_h%ns) )
+       DO  l = 0, 3
+          ALLOCATE ( surf_lsm_v(l)%pt_surface(1:surf_lsm_v(l)%ns) )
+       ENDDO 
 
 !
 !--    Allocate 2D vegetation model arrays
-       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) )
-
-       water_surface = .FALSE.
-       pave_surface  = .FALSE.
-
+!--    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_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%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_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%rad_net_l(1:surf_lsm_h%ns)        )
+       ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns)              )
+       ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)         )
+       ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns)           )
+       ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns)       )
+       ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)              )
+       ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)     )
+       ALLOCATE ( surf_lsm_h%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.
+!
+!--    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_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)%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_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)%rad_net_l(1:surf_lsm_v(l)%ns)        )
+          ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)              )
+          ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)         )
+          ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns)           )
+          ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns)       )
+          ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)              )
+          ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)     )
+          ALLOCATE ( surf_lsm_v(l)%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.
+       ENDDO
+   
 #if ! defined( __nopointer )
 !
 !--    Initial assignment of the pointers
-       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
-       t_surface => t_surface_1; t_surface_p => t_surface_2
-       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
-       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
+!--    Horizontal surfaces
+       t_soil_h    => t_soil_h_1;    t_soil_h_p    => t_soil_h_2
+       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
+!
+!--    Vertical surfaces
+       t_soil_v    => t_soil_v_1;    t_soil_v_p    => t_soil_v_2
+       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
 #endif
 
@@ -2043,5 +2487,5 @@
                                   conserve_water_content,                      &
                                   f_shortwave_incoming, field_capacity,        & 
-                                  hydraulic_conductivity,                      &
+                                  aero_resist_kray, hydraulic_conductivity,                      &
                                   lambda_surface_stable,                       &
                                   lambda_surface_unstable, leaf_area_index,    &
@@ -2087,14 +2531,16 @@
 !> temperature and water content.
 !------------------------------------------------------------------------------!
-    SUBROUTINE lsm_soil_model
+    SUBROUTINE lsm_soil_model( horizontal, l )
 
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i   !< running index
-       INTEGER(iwp) ::  j   !< running index
-       INTEGER(iwp) ::  k   !< running index
-
-       REAL(wp)     :: h_vg !< Van Genuchten coef. h
+       INTEGER(iwp) ::  k       !< running index
+       INTEGER(iwp) ::  l       !< surface-data type index indication facing
+       INTEGER(iwp) ::  m       !< running index
+
+       LOGICAL      ::  horizontal !< flag indication horizontal wall, required to set pointer accordingly
+
+       REAL(wp)     ::  h_vg !< Van Genuchten coef. h
 
        REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
@@ -2102,271 +2548,294 @@
                                                  tend           !< tendency
 
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-
-             IF ( pave_surface(j,i) )  THEN
-                rho_c_total(nzb_soil,j,i) = pave_heat_capacity
-                lambda_temp(nzb_soil)     = pave_heat_conductivity
+       TYPE(surf_type_lsm), POINTER ::  surf_m_soil
+       TYPE(surf_type_lsm), POINTER ::  surf_m_soil_p
+       TYPE(surf_type_lsm), POINTER ::  surf_t_soil
+       TYPE(surf_type_lsm), POINTER ::  surf_t_soil_p
+       TYPE(surf_type_lsm), POINTER ::  surf_tm_soil_m
+       TYPE(surf_type_lsm), POINTER ::  surf_tt_soil_m
+
+       TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable
+
+       IF ( horizontal )  THEN
+          surf           => surf_lsm_h
+
+          surf_m_soil    => m_soil_h
+          surf_m_soil_p  => m_soil_h_p
+          surf_t_soil    => t_soil_h
+          surf_t_soil_p  => t_soil_h_p
+          surf_tm_soil_m => tm_soil_h_m
+          surf_tt_soil_m => tt_soil_h_m
+       ELSE
+          surf           => surf_lsm_v(l)
+
+          surf_m_soil    => m_soil_v(l)
+          surf_m_soil_p  => m_soil_v_p(l)
+          surf_t_soil    => t_soil_v(l)
+          surf_t_soil_p  => t_soil_v_p(l)
+          surf_tm_soil_m => tm_soil_v_m(l)
+          surf_tt_soil_m => tt_soil_v_m(l)
+       ENDIF
+
+       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
+                    
+                   surf%rho_c_total(k,m) = pave_heat_capacity
+                   lambda_temp(k)        = pave_heat_conductivity   
+
+                ELSE            
+!
+!--                Calculate volumetric heat capacity of the soil, taking 
+!--                into account water content 
+                   surf%rho_c_total(k,m) = (rho_c_soil *                 &
+                                               ( 1.0_wp - surf%m_sat(m) )&
+                                               + rho_c_water * surf_m_soil%var_2d(k,m) )
+
+!
+!--                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_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) ) )
+
+                   lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +       &
+                                    lambda_h_dry
+                ENDIF
+             ENDDO
+
+!
+!--          Calculate soil heat conductivity (lambda_h) at the _stag 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)
+                ELSE
+                   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)
+
+!
+!--          Prognostic equation for soil temperature t_soil
+             tend(:) = 0.0_wp
+
+             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *&
+                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
+                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil(nzb_soil+1)            &
+                      + surf%ghf_eb(m) ) * ddz_soil_stag(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)
+
+             ENDDO
+
+             surf_t_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_t_soil%var_2d(nzb_soil:nzt_soil,m)       &
+                                               + dt_3d * ( tsc(2)              &
+                                               * tend(nzb_soil:nzt_soil)       & 
+                                               + tsc(3)                        &
+                                               * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
+
+!
+!--          Calculate t_soil tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   DO  k = nzb_soil, nzt_soil
+                      surf_tt_soil_m%var_2d(k,m) = tend(k)
+                   ENDDO
+                ELSEIF ( intermediate_timestep_count <                         &
+                         intermediate_timestep_count_max )  THEN
+                   DO  k = nzb_soil, nzt_soil
+                      surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) + 5.3125_wp        &
+                                         * surf_tt_soil_m%var_2d(k,m)
+                   ENDDO
+                ENDIF
              ENDIF
 
-             IF (  .NOT.  water_surface(j,i) )  THEN
+
+             DO  k = nzb_soil, nzt_soil
+
+!
+!--             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) )**(&
+                                                            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)&
+                                                               )**(            &
+                              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
+
+             ENDDO
+
+!
+!--          Prognostic equation for soil moisture content. Only performed,
+!--          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 
+!--             using linear interpolation. To do: replace this with
+!--             ECMWF-IFS Eq. 8.81
+                DO  k = nzb_soil, nzt_soil-1
+                    
+                   surf%lambda_w(k,m) = ( lambda_temp(k+1) +             &
+                                                lambda_temp(k) ) * 0.5_wp
+                   surf%gamma_w(k,m)  = ( gamma_temp(k+1)  +             &
+                                                gamma_temp(k) )  * 0.5_wp
+
+                ENDDO
+
+!
+!
+!--             In case of a closed bottom (= water content is conserved), 
+!--             set hydraulic conductivity to zero to that no water will be 
+!--             lost in the bottom layer.
+                IF ( conserve_water_content )  THEN
+                   surf%gamma_w(nzt_soil,m) = 0.0_wp
+                ELSE
+                   surf%gamma_w(nzt_soil,m) = gamma_temp(nzt_soil)
+                ENDIF     
+
+!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l     
+!--             * l_v)) ensures the mass conservation for water. The         
+!--             transpiration of plants equals the cumulative withdrawals by 
+!--             the roots in the soil. The scheme takes into account the 
+!--             availability of water in the soil layers as well as the root 
+!--             fraction in the respective layer. Layer with moisture below 
+!--             wilting point will not contribute, which reflects the 
+!--             preference of plants to take water from moister layers.
+!
+!--             Calculate the root extraction (ECMWF 7.69, the sum of 
+!--             root_extr = 1). The energy balance solver guarantees a 
+!--             positive transpiration, so that there is no need for an 
+!--             additional check.
+                m_total = 0.0_wp
                 DO  k = nzb_soil, nzt_soil
-
-
-                   IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
-                    
-                      rho_c_total(k,j,i) = pave_heat_capacity
-                      lambda_temp(k)     = pave_heat_conductivity   
-
-                   ELSE            
-!
-!--                   Calculate volumetric heat capacity of the soil, taking 
-!--                   into account water content 
-                      rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
-                                           + rho_c_water * m_soil(k,j,i))
-
-!
-!--                   Calculate soil heat conductivity at the center of the soil
-!--                   layers
-                      lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
-                                     lambda_h_water ** m_soil(k,j,i)
-
-                      ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i)             &
-                                                     / m_sat(j,i)))
-
-                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
-                                       lambda_h_dry
-                   ENDIF
-
+                    IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(m) )  THEN
+                       m_total = m_total + surf%root_fr(k,m) * surf_m_soil%var_2d(k,m)
+                    ENDIF
+                ENDDO  
+                 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
+                         root_extr(k) = surf%root_fr(k,m) * surf_m_soil%var_2d(k,m)  &
+                                                            / m_total
+                      ELSE
+                         root_extr(k) = 0.0_wp
+                      ENDIF
+                   ENDDO
+                ENDIF
+!
+!--             Prognostic equation for soil water content m_soil_h.
+                tend(:) = 0.0_wp
+
+                tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (       &
+                         surf_m_soil%var_2d(nzb_soil+1,m) - surf_m_soil%var_2d(nzb_soil,m) )           &
+                         * ddz_soil(nzb_soil+1) - 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)
+
+                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%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%gamma_w(k-1,m) - (root_extr(k)       &
+                             * surf%qsws_veg_eb(m) * drho_l_lv)          &
+                             ) * ddz_soil_stag(k)
                 ENDDO
-
-!
-!--             Calculate soil heat conductivity (lambda_h) at the _stag level 
-!--             using linear interpolation. For pavement surface, the
-!--             true pavement depth is considered
-                DO  k = nzb_soil, nzt_soil-1
-                   IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
-                                           .AND.  zs(k+1) > pave_depth )  THEN
-                      lambda_h(k,j,i) = ( 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)
-                   ELSE
-                      lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
-                                        * 0.5_wp
-                   ENDIF
+                tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)            &
+                                     - surf%lambda_w(nzt_soil-1,m)       &
+                                     * ( surf_m_soil%var_2d(nzt_soil,m)                    &
+                                     - surf_m_soil%var_2d(nzt_soil-1,m))                   &
+                                     * ddz_soil(nzt_soil)                      &
+                                     + 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_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_m_soil%var_2d(nzb_soil:nzt_soil,m)    &
+                                                + dt_3d * ( tsc(2) * tend(:)   &
+                                                + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
+   
+!
+!--             Account for dry soils (find a better solution here!)
+                DO  k = nzb_soil, nzt_soil
+                   IF ( surf_m_soil_p%var_2d(k,m) < 0.0_wp )  surf_m_soil_p%var_2d(k,m) = 0.0_wp
                 ENDDO
-                lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
-
-
-
-
-!
-!--             Prognostic equation for soil temperature t_soil
-                tend(:) = 0.0_wp
-
-                tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
-                          ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)  &
-                            - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
-                            + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
-
-                DO  k = nzb_soil+1, nzt_soil
-                   tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
-                             * (   lambda_h(k,j,i)                             &
-                                 * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
-                                 * ddz_soil(k+1)                               &
-                                 - lambda_h(k-1,j,i)                           &
-                                 * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
-                                 * ddz_soil(k)                                 &
-                               ) * ddz_soil_stag(k)
-
-                ENDDO
-
-                t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
-                                                  + dt_3d * ( tsc(2)           &
-                                                  * tend(nzb_soil:nzt_soil)    & 
-                                                  + tsc(3)                     &
-                                                  * tt_soil_m(:,j,i) )   
-
-!
-!--             Calculate t_soil tendencies for the next Runge-Kutta step
+
+!
+!--             Calculate m_soil tendencies for the next Runge-Kutta step
                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
                    IF ( intermediate_timestep_count == 1 )  THEN
                       DO  k = nzb_soil, nzt_soil
-                         tt_soil_m(k,j,i) = tend(k)
+                         surf_tm_soil_m%var_2d(k,m) = tend(k)
                       ENDDO
                    ELSEIF ( intermediate_timestep_count <                      &
                             intermediate_timestep_count_max )  THEN
                       DO  k = nzb_soil, nzt_soil
-                         tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
-                                         * tt_soil_m(k,j,i)
+                         surf_tm_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) + 5.3125_wp     &
+                                  * surf_tm_soil_m%var_2d(k,m)
                       ENDDO
                    ENDIF
                 ENDIF
-
-
-                DO  k = nzb_soil, nzt_soil
-
-!
-!--                Calculate soil diffusivity at the center of the soil layers
-                   lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
-                                     / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
-                                     m_wilt(j,i) ) / m_sat(j,i) )**(           &
-                                     b_ch + 2.0_wp )
-
-!
-!--                Parametrization of Van Genuchten
-                   IF ( soil_type /= 7 )  THEN
-!
-!--                   Calculate the hydraulic conductivity after Van Genuchten 
-!--                   (1980)
-                      h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
-                                 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
-                                 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
-                             )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
-
-
-                      gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
-                                      ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
-                                      1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
-                                      alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
-                                      - 1.0_wp) )**2 )                         &
-                                      / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
-                                      )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
-                                      / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
-
-!
-!--                Parametrization of Clapp & Hornberger
-                   ELSE
-                      gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
-                                      / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
-                   ENDIF
-
-                ENDDO
-
-!
-!--             Prognostic equation for soil moisture content. Only performed,
-!--             when humidity is enabled in the atmosphere and the surface type
-!--             is not pavement (implies dry soil below).
-                IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
-!
-!--                Calculate soil diffusivity (lambda_w) at the _stag level 
-!--                using linear interpolation. To do: replace this with
-!--                ECMWF-IFS Eq. 8.81
-                   DO  k = nzb_soil, nzt_soil-1
-                     
-                      lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
-                                        * 0.5_wp
-                      gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
-                                        * 0.5_wp
-
-                   ENDDO
-
-!
-!
-!--                In case of a closed bottom (= water content is conserved), 
-!--                set hydraulic conductivity to zero to that no water will be 
-!--                lost in the bottom layer.
-                   IF ( conserve_water_content )  THEN
-                      gamma_w(nzt_soil,j,i) = 0.0_wp
-                   ELSE
-                      gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
-                   ENDIF     
-
-!--                The root extraction (= root_extr * qsws_veg_eb / (rho_l     
-!--                * l_v)) ensures the mass conservation for water. The         
-!--                transpiration of plants equals the cumulative withdrawals by 
-!--                the roots in the soil. The scheme takes into account the 
-!--                availability of water in the soil layers as well as the root 
-!--                fraction in the respective layer. Layer with moisture below 
-!--                wilting point will not contribute, which reflects the 
-!--                preference of plants to take water from moister layers.
-
-!
-!--                Calculate the root extraction (ECMWF 7.69, the sum of 
-!--                root_extr = 1). The energy balance solver guarantees a 
-!--                positive transpiration, so that there is no need for an 
-!--                additional check.
-                   m_total = 0.0_wp
-                   DO  k = nzb_soil, nzt_soil
-                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
-                          m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
-                       ENDIF
-                   ENDDO  
-
-                   IF ( m_total > 0.0_wp )  THEN
-                      DO  k = nzb_soil, nzt_soil
-                         IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
-                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)      &
-                                                            / m_total
-                         ELSE
-                            root_extr(k) = 0.0_wp
-                         ENDIF
-                      ENDDO
-                   ENDIF
-
-!
-!--                Prognostic equation for soil water content m_soil.
-                   tend(:) = 0.0_wp
-
-                   tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
-                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
-                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
-                               root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
-                               + qsws_soil_eb(j,i) ) * drho_l_lv )             &
-                               * ddz_soil_stag(nzb_soil)
-
-                   DO  k = nzb_soil+1, nzt_soil-1
-                      tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
-                                - m_soil(k,j,i) ) * ddz_soil(k+1)              &
-                                - gamma_w(k,j,i)                               &
-                                - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
-                                m_soil(k-1,j,i)) * ddz_soil(k)                 &
-                                + gamma_w(k-1,j,i) - (root_extr(k)             &
-                                * qsws_veg_eb(j,i) * drho_l_lv)                &
-                                ) * ddz_soil_stag(k)
-
-                   ENDDO
-                   tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
-                                           - lambda_w(nzt_soil-1,j,i)          &
-                                           * (m_soil(nzt_soil,j,i)             &
-                                           - m_soil(nzt_soil-1,j,i))           &
-                                           * ddz_soil(nzt_soil)                &
-                                           + gamma_w(nzt_soil-1,j,i) - (       &
-                                             root_extr(nzt_soil)               &
-                                           * qsws_veg_eb(j,i) * drho_l_lv  )   &
-                                     ) * ddz_soil_stag(nzt_soil)             
-
-                   m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
-                                                   + dt_3d * ( tsc(2) * tend(:)   &
-                                                   + tsc(3) * tm_soil_m(:,j,i) )   
-   
-!
-!--                Account for dry soils (find a better solution here!)
-                   DO  k = nzb_soil, nzt_soil
-                      IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
-                   ENDDO
-
-!
-!--                Calculate m_soil tendencies for the next Runge-Kutta step
-                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                      IF ( intermediate_timestep_count == 1 )  THEN
-                         DO  k = nzb_soil, nzt_soil
-                            tm_soil_m(k,j,i) = tend(k)
-                         ENDDO
-                      ELSEIF ( intermediate_timestep_count <                   &
-                               intermediate_timestep_count_max )  THEN
-                         DO  k = nzb_soil, nzt_soil
-                            tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
-                                     * tm_soil_m(k,j,i)
-                         ENDDO
-                      ENDIF
-                   ENDIF
-
-                ENDIF
-
              ENDIF
 
-          ENDDO
+          ENDIF
+
        ENDDO
 
@@ -2386,10 +2855,19 @@
 
 #if defined( __nopointer )
-
-       t_surface    = t_surface_p
-       t_soil       = t_soil_p
+!
+!--    Horizontal surfaces
+       t_surface_h  = t_surface_h_p
+       t_soil_h     = t_soil_h_p
        IF ( humidity )  THEN
-          m_soil    = m_soil_p
-          m_liq_eb  = m_liq_eb_p
+          m_soil_h    = m_soil_h_p
+          m_liq_eb_h  = m_liq_eb_h_p
+       ENDIF
+!
+!--    Vertical surfaces
+       t_surface_v  = t_surface_v_p
+       t_soil_v     = t_soil_v_p
+       IF ( humidity )  THEN
+          m_soil_v    = m_soil_v_p
+          m_liq_eb_v  = m_liq_eb_v_p
        ENDIF
 
@@ -2399,20 +2877,39 @@
 
           CASE ( 0 )
-
-             t_surface  => t_surface_1; t_surface_p  => t_surface_2
-             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
+!
+!--          Horizontal surfaces
+             t_surface_h  => t_surface_h_1; t_surface_h_p  => t_surface_h_2
+             t_soil_h     => t_soil_h_1;    t_soil_h_p     => t_soil_h_2
              IF ( humidity )  THEN
-                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
-                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_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
              ENDIF
+!
+!--          Vertical surfaces
+             t_surface_v  => t_surface_v_1; t_surface_v_p  => t_surface_v_2
+             t_soil_v     => t_soil_v_1;    t_soil_v_p     => t_soil_v_2
+             IF ( humidity )  THEN
+                m_soil_v    => m_soil_v_1;   m_soil_v_p    => m_soil_v_2
+                m_liq_eb_v  => m_liq_eb_v_1; m_liq_eb_v_p  => m_liq_eb_v_2
+             ENDIF
+
 
 
           CASE ( 1 )
-
-             t_surface  => t_surface_2; t_surface_p  => t_surface_1
-             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
+!
+!--          Horizontal surfaces
+             t_surface_h  => t_surface_h_2; t_surface_h_p  => t_surface_h_1
+             t_soil_h     => t_soil_h_2;    t_soil_h_p     => t_soil_h_1
              IF ( humidity )  THEN
-                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
-                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
+                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
+             ENDIF
+!
+!--          Vertical surfaces
+             t_surface_v  => t_surface_v_2; t_surface_v_p  => t_surface_v_1
+             t_soil_v     => t_soil_v_2;    t_soil_v_p     => t_soil_v_1
+             IF ( humidity )  THEN
+                m_soil_v    => m_soil_v_2;   m_soil_v_p    => m_soil_v_1
+                m_liq_eb_v  => m_liq_eb_v_2; m_liq_eb_v_p  => m_liq_eb_v_1
              ENDIF
 
@@ -2445,7 +2942,8 @@
     CHARACTER (LEN=*) :: variable !< 
 
-    INTEGER(iwp) ::  i !< 
-    INTEGER(iwp) ::  j !< 
-    INTEGER(iwp) ::  k !< 
+    INTEGER(iwp) ::  i       !< 
+    INTEGER(iwp) ::  j       !< 
+    INTEGER(iwp) ::  k       !< 
+    INTEGER(iwp) ::  m       !< running index
 
     IF ( mode == 'allocate' )  THEN
@@ -2553,109 +3051,112 @@
 
           CASE ( 'c_liq*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_liq_av(j,i) = c_liq_av(j,i) + c_liq(j,i)
+             DO  m = 1, surf_lsm_h%ns
+                i   = surf_lsm_h%i(m)            
+                j   = surf_lsm_h%j(m)
+                c_liq_av(j,i) = c_liq_av(j,i) + surf_lsm_h%c_liq(m)
+             ENDDO
+
+          CASE ( 'c_soil*' )
+             DO  m = 1, surf_lsm_h%ns
+                i   = surf_lsm_h%i(m)            
+                j   = surf_lsm_h%j(m)
+                c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - surf_lsm_h%c_veg(m))
+             ENDDO
+
+          CASE ( 'c_veg*' )
+             DO  m = 1, surf_lsm_h%ns
+                i   = surf_lsm_h%i(m)            
+                j   = surf_lsm_h%j(m)
+                c_veg_av(j,i) = c_veg_av(j,i) + surf_lsm_h%c_veg(m)
+             ENDDO
+
+          CASE ( 'ghf_eb*' )
+             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)
+             ENDDO
+
+          CASE ( 'lai*' )
+             DO  m = 1, surf_lsm_h%ns
+                i   = surf_lsm_h%i(m)            
+                j   = surf_lsm_h%j(m)
+                lai_av(j,i) = lai_av(j,i) + surf_lsm_h%lai(m)
+             ENDDO
+
+          CASE ( 'm_liq_eb*' )
+             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)
+             ENDDO
+
+          CASE ( 'm_soil' )
+             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
+                   m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil_h%var_2d(k,m)
                 ENDDO
              ENDDO
 
-          CASE ( 'c_soil*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - c_veg(j,i))
-                ENDDO
-             ENDDO
-
-          CASE ( 'c_veg*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_veg_av(j,i) = c_veg_av(j,i) + c_veg(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'ghf_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   ghf_eb_av(j,i) = ghf_eb_av(j,i) + ghf_eb(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'lai*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   lai_av(j,i) = lai_av(j,i) + lai(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'm_liq_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'm_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb_soil, nzt_soil
-                      m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-
           CASE ( 'qsws_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_eb_av(j,i) = qsws_eb_av(j,i) + qsws_eb(j,i)
-                ENDDO
+             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*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) + qsws_liq_eb(j,i)
-                ENDDO
+             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*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) + qsws_soil_eb(j,i)
-                ENDDO
+             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  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) + qsws_veg_eb(j,i)
-                ENDDO
+             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)
              ENDDO
 
           CASE ( 'r_a*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   r_a_av(j,i) = r_a_av(j,i) + r_a(j,i)
-                ENDDO
+             DO  m = 1, surf_lsm_h%ns
+                i   = surf_lsm_h%i(m)            
+                j   = surf_lsm_h%j(m)
+                r_a_av(j,i) = r_a_av(j,i) + surf_lsm_h%r_a(m)
              ENDDO
 
           CASE ( 'r_s*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   r_s_av(j,i) = r_s_av(j,i) + r_s(j,i)
-                ENDDO
+             DO  m = 1, surf_lsm_h%ns
+                i   = surf_lsm_h%i(m)            
+                j   = surf_lsm_h%j(m)
+                r_s_av(j,i) = r_s_av(j,i) + surf_lsm_h%r_s(m)
              ENDDO
 
           CASE ( 'shf_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   shf_eb_av(j,i) = shf_eb_av(j,i) + shf_eb(j,i)
-                ENDDO
+             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
 
           CASE ( 't_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb_soil, nzt_soil
-                      t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil(k,j,i)
-                   ENDDO
+             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
+                   t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil_h%var_2d(k,m)
                 ENDDO
              ENDDO
@@ -2671,6 +3172,6 @@
 
           CASE ( 'c_liq*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2678,6 +3179,6 @@
 
           CASE ( 'c_soil*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2685,6 +3186,6 @@
 
           CASE ( 'c_veg*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2692,6 +3193,6 @@
 
           CASE ( 'ghf_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2699,6 +3200,6 @@
 
          CASE ( 'lai*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2706,6 +3207,6 @@
 
           CASE ( 'm_liq_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             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 )
                 ENDDO
@@ -2713,6 +3214,6 @@
 
           CASE ( 'm_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    DO  k = nzb_soil, nzt_soil
                       m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
@@ -2722,6 +3223,6 @@
 
           CASE ( 'qsws_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2729,6 +3230,6 @@
 
           CASE ( 'qsws_liq_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             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 )
                 ENDDO
@@ -2736,6 +3237,6 @@
 
           CASE ( 'qsws_soil_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             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
@@ -2743,6 +3244,6 @@
 
           CASE ( 'qsws_veg_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             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 )
                 ENDDO
@@ -2750,6 +3251,6 @@
 
           CASE ( 'r_a*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    r_a_av(j,i) = r_a_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2757,6 +3258,6 @@
 
           CASE ( 'r_s*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    r_s_av(j,i) = r_s_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
@@ -2764,6 +3265,6 @@
 
           CASE ( 't_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
                    DO  k = nzb_soil, nzt_soil
                       t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
@@ -2771,4 +3272,6 @@
                 ENDDO
              ENDDO
+!
+!-- 
 
        END SELECT
@@ -2837,8 +3340,9 @@
     CHARACTER (LEN=*) ::  variable !< 
 
-    INTEGER(iwp) ::  av !< 
-    INTEGER(iwp) ::  i  !< 
-    INTEGER(iwp) ::  j  !< 
-    INTEGER(iwp) ::  k  !< 
+    INTEGER(iwp) ::  av      !< 
+    INTEGER(iwp) ::  i       !< running index 
+    INTEGER(iwp) ::  j       !< running index
+    INTEGER(iwp) ::  k       !< running index
+    INTEGER(iwp) ::  m       !< running index
     INTEGER(iwp) ::  nzb_do  !< 
     INTEGER(iwp) ::  nzt_do  !< 
@@ -2849,15 +3353,17 @@
     REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !< 
 
+
     found = .TRUE.
 
     SELECT CASE ( TRIM( variable ) )
-
-
+!
+!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein. 
+!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged. 
        CASE ( 'c_liq*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
-                ENDDO
+             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%c_liq(m) * surf_lsm_h%c_veg(m)
              ENDDO
           ELSE
@@ -2874,8 +3380,8 @@
        CASE ( 'c_soil*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
-                ENDDO
+             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) = 1.0_wp - surf_lsm_h%c_veg(m)
              ENDDO
           ELSE
@@ -2892,8 +3398,8 @@
        CASE ( 'c_veg*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = c_veg(j,i)
-                ENDDO
+             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%c_veg(m)
              ENDDO
           ELSE
@@ -2910,8 +3416,8 @@
        CASE ( 'ghf_eb*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = ghf_eb(j,i)
-                ENDDO
+             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)
              ENDDO
           ELSE
@@ -2928,8 +3434,8 @@
        CASE ( 'lai*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = lai(j,i)
-                ENDDO
+             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%lai(m)
              ENDDO
           ELSE
@@ -2946,8 +3452,8 @@
        CASE ( 'm_liq_eb*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = m_liq_eb(j,i)
-                ENDDO
+             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)
              ENDDO
           ELSE
@@ -2964,9 +3470,9 @@
        CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO k = nzb_soil, nzt_soil
-                      local_pf(i,j,k) = m_soil(k,j,i)
-                   ENDDO
+             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
+                   local_pf(i,j,k) = m_soil_h%var_2d(k,m)
                 ENDDO
              ENDDO
@@ -2988,8 +3494,8 @@
        CASE ( 'qsws_eb*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) =  qsws_eb(j,i)
-                ENDDO
+             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)
              ENDDO
           ELSE
@@ -3006,8 +3512,8 @@
        CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
-                ENDDO
+             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)
              ENDDO
           ELSE
@@ -3024,8 +3530,8 @@
        CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
-                ENDDO
+             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)
              ENDDO
           ELSE
@@ -3042,8 +3548,8 @@
        CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
-                ENDDO
+             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
@@ -3061,8 +3567,8 @@
        CASE ( 'r_a*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = r_a(j,i)
-                ENDDO
+             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%r_a(m)
              ENDDO
           ELSE
@@ -3079,8 +3585,8 @@
        CASE ( 'r_s*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) = r_s(j,i)
-                ENDDO
+             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%r_s(m)
              ENDDO
           ELSE
@@ -3097,8 +3603,8 @@
        CASE ( 'shf_eb*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   local_pf(i,j,nzb+1) =  shf_eb(j,i)
-                ENDDO
+             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
@@ -3115,9 +3621,9 @@
        CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
           IF ( av == 0 )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO k = nzb_soil, nzt_soil
-                      local_pf(i,j,k) = t_soil(k,j,i)
-                   ENDDO
+             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
+                   local_pf(i,j,k) = t_soil_h%var_2d(k,m)
                 ENDDO
              ENDDO
@@ -3168,4 +3674,5 @@
     INTEGER(iwp) ::  j     !< 
     INTEGER(iwp) ::  k     !< 
+    INTEGER(iwp) ::  m     !< running index
 
     LOGICAL      ::  found !< 
@@ -3178,15 +3685,16 @@
 
     SELECT CASE ( TRIM( variable ) )
-
+!
+!--   Requires 3D exchange
 
       CASE ( 'm_soil' )
 
          IF ( av == 0 )  THEN
-            DO  i = nxlg, nxrg
-               DO  j = nysg, nyng
-                  DO  k = nzb_soil, nzt_soil
-                     local_pf(i,j,k) = m_soil(k,j,i)
-                  ENDDO
-               ENDDO
+            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
+                   local_pf(i,j,k) = m_soil_h%var_2d(k,m)
+                ENDDO
             ENDDO
          ELSE
@@ -3203,9 +3711,9 @@
 
          IF ( av == 0 )  THEN
-            DO  i = nxlg, nxrg
-               DO  j = nysg, nyng
-                  DO  k = nzb_soil, nzt_soil
-                     local_pf(i,j,k) = t_soil(k,j,i)
-                  ENDDO
+            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
+                  local_pf(i,j,k) = t_soil_h%var_2d(k,m)
                ENDDO
             ENDDO
@@ -3261,9 +3769,9 @@
           WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
        ENDIF
-       WRITE ( 14 )  'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb
+       WRITE ( 14 )     'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb_h
        IF ( ALLOCATED( m_liq_eb_av ) )  THEN
           WRITE ( 14 )  'm_liq_eb_av         ';  WRITE ( 14 )  m_liq_eb_av
        ENDIF
-       WRITE ( 14 )  'm_soil              ';  WRITE ( 14 )  m_soil
+       WRITE ( 14 )     'm_soil              ';  WRITE ( 14 )  m_soil_h
        IF ( ALLOCATED( m_soil_av ) )  THEN
           WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
@@ -3284,5 +3792,5 @@
           WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
        ENDIF
-       WRITE ( 14 )  't_soil              ';  WRITE ( 14 )  t_soil
+       WRITE ( 14 )     't_soil              ';  WRITE ( 14 )  t_soil_h
        IF ( ALLOCATED( t_soil_av ) )  THEN
           WRITE ( 14 )  't_soil_av           ';  WRITE ( 14 )  t_soil_av
@@ -3349,4 +3857,16 @@
           tmp_3d2   !< 
 
+    REAL(wp),                                                                  &
+       DIMENSION(1:surf_lsm_h%ns) ::                                  &
+          tmp_walltype_1d   !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ::              &
+          tmp_walltype_2d   !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ::                &
+          tmp_walltype_2d2  !< 
+
 
    IF ( initializing_actions == 'read_restart_data' )  THEN
@@ -3366,6 +3886,6 @@
             nync = nynfa(i,k) + offset_ya(i,k)
 
-
             SELECT CASE ( TRIM( field_char ) )
+
 
                 CASE ( 'c_liq_av' )
@@ -3374,5 +3894,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3382,5 +3902,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3390,5 +3910,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3398,11 +3918,11 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'm_liq_eb' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   m_liq_eb(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =        &
-                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_1d !tmp_2d
+                   m_liq_eb_h%var_1d(1:surf_lsm_h%ns)  =                       &
+                                 tmp_walltype_1d(1:surf_lsm_h%ns)
 
                 CASE ( 'lai_av' )
@@ -3411,5 +3931,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3423,8 +3943,7 @@
 
                 CASE ( 'm_soil' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
-                   m_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
-                          tmp_3d2(nzb_soil:nzt_soil,nysf-nbgp:nynf             &
-                          +nbgp,nxlf-nbgp:nxrf+nbgp)
+                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d2(:,:)
+                   m_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
+                          tmp_walltype_2d2(:,1:surf_lsm_h%ns)   
 
                 CASE ( 'm_soil_av' )
@@ -3442,5 +3961,5 @@
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                   qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3450,5 +3969,5 @@
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                   qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
                 CASE ( 'qsws_soil_eb_av' )
@@ -3473,12 +3992,11 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =       &
                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 't_soil' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   t_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
-                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,               &
-                                                nxlf-nbgp:nxrf+nbgp)
+                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d(:,:)
+                   t_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
+                         tmp_walltype_2d(:,1:surf_lsm_h%ns)      
 
                 CASE ( 't_soil_av' )
@@ -3488,5 +4006,5 @@
                    IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
                    t_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
-                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,             &
+                                    tmp_3d2(:,nysf-nbgp:nynf+nbgp,             &
                                     nxlf-nbgp:nxrf+nbgp)
 
@@ -3525,6 +4043,7 @@
        IMPLICIT NONE
 
-       INTEGER :: i  !< running index
-       INTEGER :: j  !< running index
+       INTEGER(iwp) ::  i       !< running index
+       INTEGER(iwp) ::  j       !< running index
+       INTEGER(iwp) ::  m       !< running index
 
        REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
@@ -3533,41 +4052,45 @@
 !       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
 
-
-       DO  i = nxlg, nxrg    
-          DO  j = nysg, nyng
-             IF ( water_surface(j,i) )  THEN
-
-!
-!--             Disabled: FLake parameterization. Ideally, the Charnock 
-!--             coefficient should depend on the water depth and the fetch
-!--             length
-!                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
+       DO  m = 1, surf_lsm_h%ns
+
+          i   = surf_lsm_h%i(m)            
+          j   = surf_lsm_h%j(m)
+          
+          IF ( surf_lsm_h%water_surface(m) )  THEN
+
+!
+!--          Disabled: FLake parameterization. Ideally, the Charnock 
+!--          coefficient should depend on the water depth and the fetch
+!--          length
+!             re_0 = z0(j,i) * us(j,i) / molecular_viscosity
 !        
-!                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
-!                              alpha_ch * us(j,i) / g )
-!
-!                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
-!                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
-
-!
-!--              Set minimum roughness length for u* > 0.2
-!                IF ( us(j,i) > 0.2_wp )  THEN
-!                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
-!                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
-!                ENDIF
-
-!
-!--             ECMWF IFS model parameterization after Beljaars (1994). At low
-!--             wind speed, the sea surface becomes aerodynamically smooth and
-!--             the roughness scales with the viscosity. At high wind speed, the
-!--             Charnock relation is used.
-                z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
-                          + ( alpha_ch * us(j,i)**2 / g )
-
-                z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
-                z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
-
-             ENDIF
-          ENDDO
+!             z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
+!                           alpha_ch * us(j,i) / g )
+!
+!             z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
+!             z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
+
+!
+!--           Set minimum roughness length for u* > 0.2
+!             IF ( us(j,i) > 0.2_wp )  THEN
+!                z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
+!                z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
+!             ENDIF
+
+!
+!--          ECMWF IFS model parameterization after Beljaars (1994). At low
+!--          wind speed, the sea surface becomes aerodynamically smooth and
+!--          the roughness scales with the viscosity. At high wind speed, the
+!--          Charnock relation is used.
+             surf_lsm_h%z0(m)  = ( 0.11_wp * molecular_viscosity /             &
+                                 surf_lsm_h%us(m) )  &
+                               + ( alpha_ch * surf_lsm_h%us(m)**2 / g )
+
+             surf_lsm_h%z0h(m) = 0.40_wp * molecular_viscosity /               &
+                                 surf_lsm_h%us(m)
+             surf_lsm_h%z0q(m) = 0.62_wp * molecular_viscosity /               &
+                                 surf_lsm_h%us(m)
+
+          ENDIF
        ENDDO
 
@@ -3575,54 +4098,4 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculation of specific humidity of the skin layer (surface). It is assumend
-!> that the skin is always saturated.
-!------------------------------------------------------------------------------!
-    SUBROUTINE calc_q_surface
-
-       IMPLICIT NONE
-
-       INTEGER :: i              !< running index
-       INTEGER :: j              !< running index
-       INTEGER :: k              !< running index
-
-       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
-
-       DO  i = nxlg, nxrg    
-          DO  j = nysg, nyng
-             k = nzb_s_inner(j,i)
-
-!
-!--          Calculate water vapour pressure at saturation
-             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
-                   - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
-
-!
-!--          Calculate specific humidity at saturation
-             q_s = 0.622_wp * e_s / surface_pressure
-
-             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
-
-!
-!--          Calculate specific humidity at surface
-             IF ( cloud_physics )  THEN
-                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
-                             * ( q(k+1,j,i) - ql(k+1,j,i) )
-             ELSE
-                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
-                             * q(k+1,j,i)
-             ENDIF
-
-!
-!--          Update virtual potential temperature
-             vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
-
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE calc_q_surface
-
 
  END MODULE land_surface_model_mod
