Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2214)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2232)
@@ -21,5 +21,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux; 
+! insteat, heat fluxes are directly applied in diffusion_s.
 ! 
 ! Former revisions:
@@ -103,4 +104,7 @@
 !>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
 !>
+!> 3. Temporarily large values of surface heat flux can be observed, up to
+!>    1.2 Km/s, which seem to be not realistic.
+!>
 !> @todo Check optimizations for RMA operations
 !> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
@@ -111,5 +115,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  zu, pt, pt_1, pt_2, p, ol, shf, ts, us, u, v, w, hyp, tend
+        ONLY:  zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
 
     USE cloud_parameters,                                                      &
@@ -137,5 +141,5 @@
     USE indices,                                                               &
         ONLY:  nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys,    &
-               nysg, nzb_s_inner, nzb_s_outer, nzb, nzt, nbgp
+               nysg, nzb, nzt, nbgp, wall_flags_0
 
     USE, INTRINSIC :: iso_c_binding 
@@ -157,4 +161,6 @@
     USE statistics,                                                            &
         ONLY:  hom, statistic_regions
+
+    USE surface_mod
 
                
@@ -265,4 +271,14 @@
         REAL(wp)                                   :: rtransp           !< 
     END TYPE
+!
+!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls. 
+    TYPE t_surf_vertical
+       REAL(wp), DIMENSION(:), ALLOCATABLE         :: t
+    END TYPE t_surf_vertical
+!
+!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls. 
+    TYPE t_wall_vertical
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE       :: t
+    END TYPE t_wall_vertical
 
 !-- arrays for calculation of svf and csf
@@ -320,5 +336,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface    
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface  
     
 !-- block variables needed for calculation of the plant canopy model inside the urban surface model
@@ -367,21 +383,24 @@
 !-- surface and material model variables for walls, ground, roofs
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        :: surface_types      !< array of types of wall parameters
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn                !< normalized wall layer depths (m)
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: ddz_wall           !< 1/dz_wall
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: ddz_wall_stag      !< 1/dz_wall_stag
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: dz_wall            !< wall grid spacing (center-center)
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: dz_wall_stag       !< wall grid spacing (edge-edge) 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: zw                 !< wall layer depths (m)
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf             !< wall surface temperature (K)
-    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_p           !< progn. wall surface temperature (K)
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h           !< wall surface temperature (K) at horizontal walls
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_p         !< progn. wall surface temperature (K) at horizontal walls
+
+    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v
+    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v_p
 #else
-    REAL(wp), DIMENSION(:), POINTER                :: t_surf
-    REAL(wp), DIMENSION(:), POINTER                :: t_surf_p 
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_1
-    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_2
+    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h
+    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h_p 
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_1
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_2
+
+    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v
+    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v_p
+
+    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_1
+    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_2
 #endif
     REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_av          !< average of wall surface temperature (K)
@@ -394,31 +413,22 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !-- parameters of the land, roof and wall surfaces
-    LOGICAL,  DIMENSION(:), ALLOCATABLE            :: isroof_surf        !< is the surface the part of a roof
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
 !-- parameters of the wall surfaces
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: c_surface          !< heat capacity of the wall surface skin ( J mâ2 Kâ1 )
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: lambda_surf        !< heat conductivity Î»S between air and surface ( W mâ2 Kâ1 )
-    
-!-- parameters of the walls material
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: thickness_wall     !< thickness of the wall, roof and soil layers
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: rho_c_wall         !< volumetric heat capacity of the material ( J m-3 K-1 ) (= 2.19E6)
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: lambda_h           !< heat conductivity Î»T of the material ( W m-1 K-1 )
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: roughness_wall     !< roughness relative to concrete
-    
-!-- output wall heat flux arrays
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: wshf               !< kinematic wall heat flux of sensible heat (needed for diffusion_s!<)
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: wshf_eb            !< wall heat flux of sensible heat in wall normal direction
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: wshf_eb_av         !< average of wshf_eb
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: wghf_eb            !< wall ground heat flux
-    REAL(wp), DIMENSION(:), ALLOCATABLE            :: wghf_eb_av         !< average of wghf_eb
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall             !< Wall temperature (K)
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall_av          !< Average of t_wall
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall_p           !< Prog. wall temperature (K)
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_wall_h             !< Wall temperature (K)
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_wall_h_av          !< Average of t_wall
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_wall_h_p           !< Prog. wall temperature (K)
+
+    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v             !< Wall temperature (K)
+    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av          !< Average of t_wall
+    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_p           !< Prog. wall temperature (K)
 #else
-    REAL(wp), DIMENSION(:,:), POINTER              :: t_wall, t_wall_p
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall_av, t_wall_1, t_wall_2
+    REAL(wp), DIMENSION(:,:), POINTER                :: t_wall_h, t_wall_h_p
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_av, t_wall_h_1, t_wall_h_2
+
+    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_wall_v, t_wall_v_p
+    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av, t_wall_v_1, t_wall_v_2
 #endif
 
@@ -488,10 +498,5 @@
        MODULE PROCEDURE usm_swap_timelevel
     END INTERFACE usm_swap_timelevel
-    
-    INTERFACE usm_wall_heat_flux
-       MODULE PROCEDURE usm_wall_heat_flux
-       MODULE PROCEDURE usm_wall_heat_flux_ij
-    END INTERFACE usm_wall_heat_flux
-    
+        
     INTERFACE usm_write_restart_data
        MODULE PROCEDURE usm_write_restart_data
@@ -508,5 +513,4 @@
            usm_energy_balance_land, usm_energy_balance_wall, nrefsteps,        &
            usm_init_urban_surface, usm_radiation, usm_read_restart_data,       &
-           usm_wall_heat_flux,                                                 &
            usm_surface_energy_balance, usm_material_heat_model,                &
            usm_swap_timelevel, usm_check_data_output, usm_average_3d_data,     &
@@ -531,7 +535,9 @@
         IMPLICIT NONE
        
-        INTEGER(iwp)                            :: i, j, k, d, l, ir, jr, ids
-        INTEGER(iwp)                            :: nzubl, nzutl, isurf, ipcgb
-        INTEGER(iwp)                            :: procid
+        INTEGER(iwp) :: i, j, k, d, l, ir, jr, ids, m
+        INTEGER(iwp) :: k_topo     !< vertical index indicating topography top for given (j,i)
+        INTEGER(iwp) :: k_topo2    !< vertical index indicating topography top for given (j,i)
+        INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb
+        INTEGER(iwp) :: procid
 
         
@@ -543,8 +549,24 @@
         CALL location_message( '', .TRUE. )
         CALL location_message( '    allocation of needed arrays', .TRUE. )
-!--     find nzub, nzut, nzu
-        nzubl = minval(nzb_s_inner(nys:nyn,nxl:nxr))
-        nzutl = maxval(nzb_s_inner(nys:nyn,nxl:nxr))
+!
+!--     Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be 
+!--     removed later). The following contruct finds the lowest / largest index
+!--     for any upward-facing wall (see bit 12). 
+        nzubl = MINVAL(                                                        &
+                    MAXLOC(                                                    &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
+                               ), DIM = 1                                      &
+                          ) - 1                                                & 
+                            )
+        nzutl = MAXVAL(                                                        &
+                   MAXLOC(                                                     &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
+                               ), DIM = 1                                      &
+                          ) - 1                                                &
+                            )
         nzubl = max(nzubl,nzb)
+
         
         IF ( plant_canopy )  THEN
@@ -559,8 +581,14 @@
             DO i = nxl, nxr
                 DO j = nys, nyn
+!
+!--                 Find topography top index
+                    k_topo = MAXLOC( MERGE( 1, 0,                              &
+                                             BTEST( wall_flags_0(:,j,i), 12 )  &
+                                           ), DIM = 1                          &
+                                   ) - 1
                     DO k = nzt+1, 0, -1
                         IF ( lad_s(k,j,i) /= 0.0_wp )  THEN
 !--                         we are at the top of the pcs
-                            pct(j,i) = k + nzb_s_inner(j,i)
+                            pct(j,i) = k + k_topo
                             pch(j,i) = k
                             npcbl = npcbl + pch(j,i)
@@ -601,27 +629,23 @@
         CALL location_message( '    calculation of indices for surfaces', .TRUE. )
         nsurfl = 0
-!--     calculate land surface and roof
-        startland = nsurfl+1
-        nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1)
-        endland = nsurfl
-        nlands = endland-startland+1
-
-!--     calculation of the walls
+!
+!--     Number of land- and roof surfaces. Note, since horizontal surface elements
+!--     are already counted in surface_mod, in case be simply reused here.
+        startland = 1
+        nsurfl    = surf_usm_h%ns
+        endland   = nsurfl
+        nlands    = endland-startland+1
+
+!
+!--     Number of vertical surfaces. As vertical surfaces are already 
+!--     counted in surface mod, it can be reused here. 
         startwall = nsurfl+1
-        DO i = nxl, nxr
-            DO j = nys, nyn
-!--             test for walls
-!--             (we don't use array flags because it isn't calculated in case of masking_method=.T.)
-                DO ids = 1, 4  !-- four wall directions
-                    jr = min(max(j-jdir(ids),0),ny)
-                    ir = min(max(i-idir(ids),0),nx)
-                    nsurfl = nsurfl + max(0, nzb_s_inner(jr,ir)-nzb_s_inner(j,i))
-                ENDDO
-            ENDDO
-        ENDDO
+        nsurfl = nsurfl + surf_usm_v(0)%ns + surf_usm_v(1)%ns +        &
+                          surf_usm_v(2)%ns + surf_usm_v(3)%ns 
         endwall = nsurfl
         nwalls = endwall-startwall+1
-        
-!--     range of energy balance surfaces
+
+        
+!--     range of energy balance surfaces  ! will be treated separately by surf_usm_h and surf_usm_v
         nenergy = 0
         IF ( usm_energy_balance_land )  THEN
@@ -641,5 +665,5 @@
 !--     block of virtual surfaces
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!--     calculate sky surfaces
+!--     calculate sky surfaces  ! not used so far!
         startsky = nsurfl+1
         nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1)
@@ -656,14 +680,24 @@
         ijdb = RESHAPE( (/ nxl,nxr,nyn,nyn,nxl,nxr,nys,nys,nxr,nxr,nys,nyn,nxl,nxl,nys,nyn /), (/4, 4/) )
 !--     calulation of the free borders of the domain
-        DO ids = 6,9
-            IF ( isborder(ids) )  THEN
-!--             free border of the domain in direction ids
-                DO i = ijdb(1,ids), ijdb(2,ids)
-                    DO j = ijdb(3,ids), ijdb(4,ids)
-                        k = nzut - max(nzb_s_inner(j,i), nzb_s_inner(j-jdir(ids),i-idir(ids)))
-                        nsurfl = nsurfl + k
-                    ENDDO
-                ENDDO
-            ENDIF
+        DO  ids = 6,9
+           IF ( isborder(ids) )  THEN
+!--           free border of the domain in direction ids
+              DO  i = ijdb(1,ids), ijdb(2,ids)
+                 DO  j = ijdb(3,ids), ijdb(4,ids)
+
+                    k_topo =   MAXLOC( MERGE( 1, 0,                            &
+                                             BTEST( wall_flags_0(:,j,i), 12 )  &
+                                            ), DIM = 1                         &
+                                     ) - 1
+                    k_topo2 =  MAXLOC( MERGE( 1, 0,                            &
+                                             BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 )  &
+                                            ), DIM = 1                         &
+                                     ) - 1 
+
+                    k = nzut - MAX( k_topo, k_topo2 )
+                    nsurfl = nsurfl + k
+                 ENDDO
+              ENDDO
+           ENDIF
         ENDDO
         
@@ -676,5 +710,11 @@
             DO i = nxl, nxr
                 DO j = nys, nyn
-                    DO k = nzb_s_inner(j,i)+1, pct(j,i)
+!
+!--                 Find topography top index
+                    k_topo = MAXLOC( MERGE( 1, 0,                              &
+                                             BTEST( wall_flags_0(:,j,i), 12 )  &
+                                          ), DIM = 1                           &
+                                   ) - 1
+                    DO k = k_topo + 1, pct(j,i)
                         ipcgb = ipcgb + 1
                         gridpcbl(k,j,i) = ipcgb
@@ -689,5 +729,5 @@
 
 !--     fill surfl
-        ALLOCATE(surfl(4,nsurfl))
+        ALLOCATE(surfl(5,nsurfl))
         isurf = 0
         
@@ -695,7 +735,10 @@
         DO i = nxl, nxr
             DO j = nys, nyn
-                isurf = isurf + 1
-                k = nzb_s_inner(j,i)+1
-                surfl(:,isurf) = (/iroof,k,j,i/)
+               DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+                  k = surf_usm_h%k(m)
+
+                  isurf = isurf + 1
+                  surfl(:,isurf) = (/iroof,k,j,i,m/)
+               ENDDO
             ENDDO
         ENDDO
@@ -704,12 +747,32 @@
         DO i = nxl, nxr
             DO j = nys, nyn
-                DO ids = 1, 4  !> four wall directions
-                    jr = min(max(j-jdir(ids),0),ny)
-                    ir = min(max(i-idir(ids),0),nx)
-                    DO k = nzb_s_inner(j,i)+1, nzb_s_inner(jr,ir)
-                        isurf = isurf + 1
-                        surfl(:,isurf) = (/ids,k,j,i/)
-                    ENDDO
-                ENDDO
+               l = 0
+               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
+                  k = surf_usm_v(l)%k(m)
+
+                  isurf          = isurf + 1
+                  surfl(:,isurf) = (/2,k,j,i,m/)
+               ENDDO
+               l = 1
+               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
+                  k = surf_usm_v(l)%k(m)
+
+                  isurf          = isurf + 1
+                  surfl(:,isurf) = (/1,k,j,i,m/)
+               ENDDO
+               l = 2
+               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
+                  k = surf_usm_v(l)%k(m)
+
+                  isurf          = isurf + 1
+                  surfl(:,isurf) = (/4,k,j,i,m/)
+               ENDDO
+               l = 3
+               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
+                  k = surf_usm_v(l)%k(m)
+
+                  isurf          = isurf + 1
+                  surfl(:,isurf) = (/3,k,j,i,m/)
+               ENDDO
             ENDDO
         ENDDO
@@ -720,5 +783,5 @@
                 isurf = isurf + 1
                 k = nzut
-                surfl(:,isurf) = (/isky,k,j,i/)
+                surfl(:,isurf) = (/isky,k,j,i,-1/)
             ENDDO
         ENDDO
@@ -730,7 +793,15 @@
                 DO i = ijdb(1,ids), ijdb(2,ids)
                     DO j = ijdb(3,ids), ijdb(4,ids)
-                        DO k = max(nzb_s_inner(j,i),nzb_s_inner(j-jdir(ids),i-idir(ids)))+1, nzut
+                        k_topo =   MAXLOC( MERGE( 1, 0,                        &
+                                             BTEST( wall_flags_0(:,j,i), 12 )  &
+                                                ), DIM = 1                     &
+                                         ) - 1
+                        k_topo2 =  MAXLOC( MERGE( 1, 0,                        &
+                                             BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 )  &
+                                                ), DIM = 1                     &
+                                         ) - 1 
+                        DO k = MAX(k_topo,k_topo2)+1, nzut
                             isurf = isurf + 1
-                            surfl(:,isurf) = (/ids,k,j,i/)
+                            surfl(:,isurf) = (/ids,k,j,i,-1/)
                         ENDDO
                     ENDDO
@@ -755,8 +826,8 @@
         surfstart(numprocs) = k
         nsurf = k
-        ALLOCATE(surf(4,nsurf))
+        ALLOCATE(surf(5,nsurf))
         
 #if defined( __parallel )
-        CALL MPI_AllGatherv(surfl, nsurfl*4, MPI_INTEGER, surf, nsurfs*4, surfstart*4, MPI_INTEGER, comm2d, ierr)
+        CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, surfstart*5, MPI_INTEGER, comm2d, ierr)
 #else
         surf = surfl
@@ -787,6 +858,28 @@
         ALLOCATE( surfouts(nsurf) ) !TODO: global surfaces without virtual
         ALLOCATE( surfoutl(nsurf) ) !TODO: global surfaces without virtual
-        ALLOCATE( surfhf(startenergy:endenergy) )
-        ALLOCATE( rad_net_l(startenergy:endenergy) )
+
+
+
+!
+!--     Allocate radiation arrays which are part of the new data type. 
+!--     For horizontal surfaces.
+        ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns)    )
+        ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
+!
+!--  New
+        ALLOCATE( surf_usm_h%rad_in_sw(1:surf_usm_h%ns)  )
+        ALLOCATE( surf_usm_h%rad_out_sw(1:surf_usm_h%ns) )
+        ALLOCATE( surf_usm_h%rad_in_lw(1:surf_usm_h%ns)  )
+        ALLOCATE( surf_usm_h%rad_out_lw(1:surf_usm_h%ns) )
+!
+!--     For vertical surfaces
+        DO  l = 0, 3
+           ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns)    )
+           ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
+           ALLOCATE( surf_usm_v(l)%rad_in_sw(1:surf_usm_v(l)%ns)  )
+           ALLOCATE( surf_usm_v(l)%rad_out_sw(1:surf_usm_v(l)%ns) )
+           ALLOCATE( surf_usm_v(l)%rad_in_lw(1:surf_usm_v(l)%ns)  )
+           ALLOCATE( surf_usm_v(l)%rad_out_lw(1:surf_usm_v(l)%ns) )
+        ENDDO
 
 !--     Wall surface model
@@ -794,5 +887,8 @@
        
 !--     allocate array of wall types and wall parameters
-        ALLOCATE ( surface_types(startenergy:endenergy) )
+        ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
+        DO  l = 0, 3
+           ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
+        ENDDO
        
 !--     broadband albedo of the land, roof and wall surface
@@ -802,58 +898,131 @@
         ALLOCATE ( albedo_surf(nsurfl) )
         albedo_surf = 1.0_wp
+        ALLOCATE ( surf_usm_h%albedo_surf(1:surf_usm_h%ns) )
+        DO  l = 0, 3
+           ALLOCATE( surf_usm_v(l)%albedo_surf(1:surf_usm_v(l)%ns) )
+        ENDDO
        
-!--     wall and roof surface parameters
-        ALLOCATE ( isroof_surf(startenergy:endenergy) )
+!--     wall and roof surface parameters. First for horizontal surfaces
         ALLOCATE ( emiss_surf(startenergy:endenergy) )
-        ALLOCATE ( lambda_surf(startenergy:endenergy) )
-        ALLOCATE ( c_surface(startenergy:endenergy) )
-        ALLOCATE ( roughness_wall(startenergy:endenergy) )
+
+        ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns)    )
+        ALLOCATE ( surf_usm_h%emiss_surf(1:surf_usm_h%ns)     )
+        ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns)    )
+        ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns)      )
+        ALLOCATE ( surf_usm_h%roughness_wall(1:surf_usm_h%ns) )
+!
+!--     For vertical surfaces.
+        DO  l = 0, 3
+           ALLOCATE ( surf_usm_v(l)%emiss_surf(1:surf_usm_v(l)%ns)     )
+           ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns)    )
+           ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns)      )
+           ALLOCATE ( surf_usm_v(l)%roughness_wall(1:surf_usm_v(l)%ns) )
+        ENDDO
        
-!--     allocate wall and roof material parameters
-        ALLOCATE ( thickness_wall(startenergy:endenergy) )
-        ALLOCATE ( lambda_h(nzb_wall:nzt_wall,startenergy:endenergy) )
-        ALLOCATE ( rho_c_wall(nzb_wall:nzt_wall,startenergy:endenergy) )
-
-!--     allocate wall and roof layers sizes
+!--     allocate wall and roof material parameters. First for horizontal surfaces
+        ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns)               )
+        ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns)   )
+        ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
+!
+!--     For vertical surfaces.
+        DO  l = 0, 3
+           ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns)               )
+           ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)   )
+           ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
+        ENDDO
+
+!--     allocate wall and roof layers sizes. For horizontal surfaces.
         ALLOCATE ( zwn(nzb_wall:nzt_wall) )
-        ALLOCATE ( dz_wall(nzb_wall:nzt_wall+1, startenergy:endenergy) )
-        ALLOCATE ( ddz_wall(nzb_wall:nzt_wall+1, startenergy:endenergy) )
-        ALLOCATE ( dz_wall_stag(nzb_wall:nzt_wall, startenergy:endenergy) )
-        ALLOCATE ( ddz_wall_stag(nzb_wall:nzt_wall, startenergy:endenergy) )
-        ALLOCATE ( zw(nzb_wall:nzt_wall, startenergy:endenergy) )
-
-!--     allocate wall and roof temperature arrays
+        ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
+        ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)    )
+        ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
+        ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
+        ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns)            )
+!
+!--     For vertical surfaces.
+        DO  l = 0, 3
+           ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
+           ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)    )
+           ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
+           ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
+           ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)            )
+        ENDDO
+
+!--     allocate wall and roof temperature arrays, for horizontal walls
 #if defined( __nopointer )
-        ALLOCATE ( t_surf(startenergy:endenergy) )
-        ALLOCATE ( t_surf_p(startenergy:endenergy) )
-        ALLOCATE ( t_wall(nzb_wall:nzt_wall+1,startenergy:endenergy) )
-        ALLOCATE ( t_wall_p(nzb_wall:nzt_wall+1,startenergy:endenergy) )
+        ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
+        ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
+        ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
+        ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
+
+        ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
+        ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
+        ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
+        ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
 #else
-        ALLOCATE ( t_surf_1(startenergy:endenergy) )
-        ALLOCATE ( t_surf_2(startenergy:endenergy) )
-        ALLOCATE ( t_wall_1(nzb_wall:nzt_wall+1,startenergy:endenergy) )
-        ALLOCATE ( t_wall_2(nzb_wall:nzt_wall+1,startenergy:endenergy) )
+        ALLOCATE ( t_surf_h_1(1:surf_usm_h%ns) )
+        ALLOCATE ( t_surf_h_2(1:surf_usm_h%ns) )
+        ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
+        ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
 
 !--     initial assignment of the pointers
-        t_wall    => t_wall_1;    t_wall_p    => t_wall_2
-        t_surf => t_surf_1; t_surf_p => t_surf_2
+        t_wall_h    => t_wall_h_1;    t_wall_h_p    => t_wall_h_2
+        t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2
 #endif
 
-!--     allocate intermediate timestep arrays
-        ALLOCATE ( tt_surface_m(startenergy:endenergy) )
-        ALLOCATE ( tt_wall_m(nzb_wall:nzt_wall+1,startenergy:endenergy) )
-
-!--     allocate wall heat flux output array
-        ALLOCATE ( wshf(startwall:endwall) )
-        ALLOCATE ( wshf_eb(startenergy:endenergy) )
-        ALLOCATE ( wghf_eb(startenergy:endenergy) )
-
-!--     set inital values for prognostic quantities
-        tt_surface_m = 0.0_wp
-        tt_wall_m    = 0.0_wp
-
-        wshf = 0.0_wp
-        wshf_eb = 0.0_wp
-        wghf_eb = 0.0_wp
+!--     allocate wall and roof temperature arrays, for vertical walls
+#if defined( __nopointer )
+        DO  l = 0, 3
+           ALLOCATE ( t_surf_v(l)%t(1:surf_usm_v(l)%ns) )
+           ALLOCATE ( t_surf_v(l)%t_p(1:surf_usm_v(l)%ns) )
+           ALLOCATE ( t_wall_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
+           ALLOCATE ( t_wall_v(l)%t_p(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
+        ENDDO
+#else
+        DO  l = 0, 3
+           ALLOCATE ( t_surf_v_1(l)%t(1:surf_usm_v(l)%ns) )
+           ALLOCATE ( t_surf_v_2(l)%t(1:surf_usm_v(l)%ns) )
+           ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
+           ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
+        ENDDO
+!
+!--     initial assignment of the pointers
+        t_wall_v    => t_wall_v_1;    t_wall_v_p    => t_wall_v_2
+        t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2
+#endif
+!
+!--     Allocate intermediate timestep arrays. For horizontal surfaces.
+        ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns)                  )
+        ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
+!
+!--     Set inital values for prognostic quantities
+        IF ( ALLOCATED( surf_usm_h%tt_surface_m ) )  surf_usm_h%tt_surface_m = 0.0_wp
+        IF ( ALLOCATED( surf_usm_h%tt_wall_m    ) )  surf_usm_h%tt_wall_m    = 0.0_wp
+!
+!--     Now, for vertical surfaces
+        DO  l = 0, 3
+           ALLOCATE ( surf_usm_v(l)%tt_surface_m(1:surf_usm_v(l)%ns)                  )
+           ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
+           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) )  surf_usm_v(l)%tt_surface_m = 0.0_wp
+           IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m    ) )  surf_usm_v(l)%tt_wall_m    = 0.0_wp
+        ENDDO
+
+!--     allocate wall heat flux output array and set initial values. For horizontal surfaces
+!         ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns)    )  !can be removed 
+        ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
+        ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
+        IF ( ALLOCATED( surf_usm_h%wshf    ) )  surf_usm_h%wshf    = 0.0_wp
+        IF ( ALLOCATED( surf_usm_h%wshf_eb ) )  surf_usm_h%wshf_eb = 0.0_wp
+        IF ( ALLOCATED( surf_usm_h%wghf_eb ) )  surf_usm_h%wghf_eb = 0.0_wp
+!
+!--     Now, for vertical surfaces
+        DO  l = 0, 3
+!            ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns)    )    ! can be removed
+           ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
+           ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
+           IF ( ALLOCATED( surf_usm_v(l)%wshf    ) )  surf_usm_v(l)%wshf    = 0.0_wp
+           IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) )  surf_usm_v(l)%wshf_eb = 0.0_wp
+           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) )  surf_usm_v(l)%wghf_eb = 0.0_wp
+        ENDDO
         
     END SUBROUTINE usm_allocate_urban_surface
@@ -874,5 +1043,5 @@
         CHARACTER (len=*), INTENT(IN) :: variable
   
-        INTEGER(iwp)                                       :: i, j, k, l, ids, iwl,istat
+        INTEGER(iwp)                                       :: i, j, k, l, m, ids, iwl,istat
         CHARACTER (len=varnamelength)                      :: var, surfid
         INTEGER(iwp), PARAMETER                            :: nd = 5
@@ -910,22 +1079,40 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    IF ( .NOT.  ALLOCATED(rad_net_av) )  THEN
-                        ALLOCATE( rad_net_av(startenergy:endenergy) )
-                        rad_net_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
+                        ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
+                        surf_usm_h%rad_net_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%rad_net_av = 0.0_wp
+                       ENDIF
+                    ENDDO
                     
                 CASE ( 'usm_rad_insw' )
 !--                 array of sw radiation falling to surface after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
-                        ALLOCATE( surfinsw_av(startenergy:endenergy) )
-                        surfinsw_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfinsw_av) )  THEN
+                        ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) )
+                        surf_usm_h%surfinsw_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinsw_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%surfinsw_av = 0.0_wp
+                       ENDIF
+                    ENDDO
                                     
                 CASE ( 'usm_rad_inlw' )
 !--                 array of lw radiation falling to surface after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
-                        ALLOCATE( surfinlw_av(startenergy:endenergy) )
-                        surfinlw_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfinlw_av) )  THEN
+                        ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) )
+                        surf_usm_h%surfinlw_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinlw_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%surfinlw_av = 0.0_wp
+                       ENDIF
+                    ENDDO
 
                 CASE ( 'usm_rad_inswdir' )
@@ -952,5 +1139,5 @@
                 CASE ( 'usm_rad_inlwdif' )
 !--                 array of sw radiation falling to surface after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
+                   IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
                         ALLOCATE( surfinlwdif_av(startenergy:endenergy) )
                         surfinlwdif_av = 0.0_wp
@@ -977,5 +1164,4 @@
                         surfoutlw_av = 0.0_wp
                     ENDIF
-
                 CASE ( 'usm_rad_ressw' )
 !--                 array of residua of sw radiation absorbed in surface after last reflection
@@ -984,5 +1170,5 @@
                         surfins_av = 0.0_wp
                     ENDIF
-                                    
+                                   
                 CASE ( 'usm_rad_reslw' )
 !--                 array of residua of lw radiation absorbed in surface after last reflection
@@ -994,37 +1180,67 @@
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfhf_av) )  THEN
-                        ALLOCATE( surfhf_av(startenergy:endenergy) )
-                        surfhf_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
+                        ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
+                        surf_usm_h%surfhf_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%surfhf_av = 0.0_wp
+                       ENDIF
+                    ENDDO
 
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces
 !--                 land surfaces
-                    IF ( .NOT.  ALLOCATED(wshf_eb_av) )  THEN
-                        ALLOCATE( wshf_eb_av(startenergy:endenergy) )
-                        wshf_eb_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
+                        ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
+                        surf_usm_h%wshf_eb_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%wshf_eb_av = 0.0_wp
+                       ENDIF
+                    ENDDO
 
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(wghf_eb_av) )  THEN
-                        ALLOCATE( wghf_eb_av(startenergy:endenergy) )
-                        wghf_eb_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
+                        ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
+                        surf_usm_h%wghf_eb_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%wghf_eb_av = 0.0_wp
+                       ENDIF
+                    ENDDO
 
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    IF ( .NOT.  ALLOCATED(t_surf_av) )  THEN
-                        ALLOCATE( t_surf_av(startenergy:endenergy) )
-                        t_surf_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
+                        ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
+                        surf_usm_h%t_surf_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%t_surf_av = 0.0_wp
+                       ENDIF
+                    ENDDO
 
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for iwl layer of walls and land
-                    IF ( .NOT.  ALLOCATED(t_wall_av) )  THEN
-                        ALLOCATE( t_wall_av(nzb_wall:nzt_wall,startenergy:endenergy) )
-                        t_wall_av = 0.0_wp
+                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
+                        ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
+                        surf_usm_h%t_wall_av = 0.0_wp
                     ENDIF
+                    DO  l = 0, 3
+                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
+                           ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
+                           surf_usm_v(l)%t_wall_av = 0.0_wp
+                       ENDIF
+                    ENDDO
 
                CASE DEFAULT
@@ -1039,8 +1255,15 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            rad_net_av(l) = rad_net_av(l) + rad_net_l(l)
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%rad_net_av(m) =                              &
+                                          surf_usm_h%rad_net_av(m) +           &
+                                          surf_usm_h%rad_net_l(m)
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%rad_net_av(m) =                        &
+                                          surf_usm_v(l)%rad_net_av(m) +        &
+                                          surf_usm_v(l)%rad_net_l(m)
+                       ENDDO
                     ENDDO
                     
@@ -1085,4 +1308,5 @@
                         ENDIF
                     ENDDO
+
                     
                 CASE ( 'usm_rad_inlwdif' )
@@ -1090,10 +1314,19 @@
                     DO l = startenergy, endenergy
                         IF ( surfl(id,l) == ids )  THEN
+                            surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
+                                                surfinswdir(l) - surfinswdif(l)
+                        ENDIF
+                    ENDDO
+!                     
+                CASE ( 'usm_rad_inlwref' )
+!--                 array of lw radiation falling to surface from reflections
+                    DO l = startenergy, endenergy
+                        IF ( surfl(id,l) == ids )  THEN
                             surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
                         ENDIF
                     ENDDO
                     
-                CASE ( 'usm_rad_inlwref' )
-!--                 array of lw radiation falling to surface from reflections
+                CASE ( 'usm_rad_outsw' )
+!--                 array of sw radiation emitted from surface after i-th reflection
                     DO l = startenergy, endenergy
                         IF ( surfl(id,l) == ids )  THEN
@@ -1103,17 +1336,9 @@
                     ENDDO
                     
-                CASE ( 'usm_rad_outsw' )
-!--                 array of sw radiation emitted from surface after i-th reflection
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
-                        ENDIF
-                    ENDDO
-                    
                 CASE ( 'usm_rad_outlw' )
 !--                 array of lw radiation emitted from surface after i-th reflection
                     DO l = startenergy, endenergy
                         IF ( surfl(id,l) == ids )  THEN
-                            surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
+                            surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
                         ENDIF
                     ENDDO
@@ -1123,5 +1348,5 @@
                     DO l = startenergy, endenergy
                         IF ( surfl(id,l) == ids )  THEN
-                            surfins_av(l) = surfins_av(l) + surfins(l)
+                            surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
                         ENDIF
                     ENDDO
@@ -1131,5 +1356,5 @@
                     DO l = startenergy, endenergy
                         IF ( surfl(id,l) == ids )  THEN
-                            surfinl_av(l) = surfinl_av(l) + surfinl(l)
+                            surfins_av(l) = surfins_av(l) + surfins(l)
                         ENDIF
                     ENDDO
@@ -1137,40 +1362,75 @@
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            surfhf_av(l) = surfhf_av(l) + surfhf(l)
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%surfhf_av(m) =                               &
+                                          surf_usm_h%surfhf_av(m) +            &
+                                          surf_usm_h%surfhf(m)
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%surfhf_av(m) =                         &
+                                          surf_usm_v(l)%surfhf_av(m) +         &
+                                          surf_usm_v(l)%surfhf(m)
+                       ENDDO
                     ENDDO
                     
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            wshf_eb_av(l) = wshf_eb_av(l) + wshf_eb(l)
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%wshf_eb_av(m) =                              &
+                                          surf_usm_h%wshf_eb_av(m) +           &
+                                          surf_usm_h%wshf_eb(m)
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%wshf_eb_av(m) =                        &
+                                          surf_usm_v(l)%wshf_eb_av(m) +        &
+                                          surf_usm_v(l)%wshf_eb(m)
+                       ENDDO
                     ENDDO
                     
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            wghf_eb_av(l) = wghf_eb_av(l) + wghf_eb(l)
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%wghf_eb_av(m) =                              &
+                                          surf_usm_h%wghf_eb_av(m) +           &
+                                          surf_usm_h%wghf_eb(m)
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%wghf_eb_av(m) =                        &
+                                          surf_usm_v(l)%wghf_eb_av(m) +        &
+                                          surf_usm_v(l)%wghf_eb(m)
+                       ENDDO
                     ENDDO
                     
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            t_surf_av(l) = t_surf_av(l) + t_surf(l)
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%t_surf_av(m) =                               & 
+                                          surf_usm_h%t_surf_av(m) +            &
+                                          t_surf_h(m)
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%t_surf_av(m) =                         &
+                                          surf_usm_v(l)%t_surf_av(m) +         &
+                                          t_surf_v(l)%t(m)
+                       ENDDO
                     ENDDO
                     
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for  iwl layer of walls and land
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            t_wall_av(iwl, l) = t_wall_av(iwl,l) + t_wall(iwl, l)
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%t_wall_av(iwl,m) =                           &
+                                          surf_usm_h%t_wall_av(iwl,m) +        &
+                                          t_wall_h(iwl,m)
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
+                                          surf_usm_v(l)%t_wall_av(iwl,m) +     &
+                                          t_wall_v(l)%t(iwl,m)
+                       ENDDO
                     ENDDO
                     
@@ -1186,8 +1446,15 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            rad_net_av(l) = rad_net_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%rad_net_av(m) =                              &
+                                          surf_usm_h%rad_net_av(m) /           &
+                                          REAL( average_count_3d, kind=wp )
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%rad_net_av(m) =                        &
+                                          surf_usm_v(l)%rad_net_av(m) /        &
+                                          REAL( average_count_3d, kind=wp )
+                       ENDDO
                     ENDDO
                     
@@ -1199,5 +1466,5 @@
                         ENDIF
                     ENDDO
-                                    
+                             
                 CASE ( 'usm_rad_inlw' )
 !--                 array of lw radiation falling to surface after i-th reflection
@@ -1207,5 +1474,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_inswdir' )
 !--                 array of direct sw radiation falling to surface from sun
@@ -1215,5 +1482,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_inswdif' )
 !--                 array of difusion sw radiation falling to surface from sky and borders of the domain
@@ -1223,5 +1490,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_inswref' )
 !--                 array of sw radiation falling to surface from reflections
@@ -1231,5 +1498,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_inlwdif' )
 !--                 array of sw radiation falling to surface after i-th reflection
@@ -1239,5 +1506,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_inlwref' )
 !--                 array of lw radiation falling to surface from reflections
@@ -1247,5 +1514,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_outsw' )
 !--                 array of sw radiation emitted from surface after i-th reflection
@@ -1255,5 +1522,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_outlw' )
 !--                 array of lw radiation emitted from surface after i-th reflection
@@ -1263,5 +1530,5 @@
                         ENDIF
                     ENDDO
-
+                    
                 CASE ( 'usm_rad_ressw' )
 !--                 array of residua of sw radiation absorbed in surface after last reflection
@@ -1282,40 +1549,75 @@
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            surfhf_av(l) = surfhf_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%surfhf_av(m) =                               &
+                                          surf_usm_h%surfhf_av(m) /            &
+                                          REAL( average_count_3d, kind=wp )
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%surfhf_av(m) =                         &
+                                          surf_usm_v(l)%surfhf_av(m) /         &
+                                          REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ENDDO
+                    
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            wshf_eb_av(l) = wshf_eb_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%wshf_eb_av(m) =                              &
+                                          surf_usm_h%wshf_eb_av(m) /           &
+                                          REAL( average_count_3d, kind=wp )
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%wshf_eb_av(m) =                        &
+                                          surf_usm_v(l)%wshf_eb_av(m) /        &
+                                          REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ENDDO
+                    
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            wghf_eb_av(l) = wghf_eb_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%wghf_eb_av(m) =                              &
+                                          surf_usm_h%wghf_eb_av(m) /           &
+                                          REAL( average_count_3d, kind=wp )
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%wghf_eb_av(m) =                        &
+                                          surf_usm_v(l)%wghf_eb_av(m) /        &
+                                          REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ENDDO
+                    
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            t_surf_av(l) = t_surf_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%t_surf_av(m) =                               & 
+                                          surf_usm_h%t_surf_av(m) /            &
+                                          REAL( average_count_3d, kind=wp )
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%t_surf_av(m) =                         &
+                                          surf_usm_v(l)%t_surf_av(m) /         &
+                                          REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ENDDO
+                    
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for  iwl layer of walls and land
-                    DO l = startenergy, endenergy
-                        IF ( surfl(id,l) == ids )  THEN
-                            t_wall_av(iwl, l) = t_wall_av(iwl,l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
+                    DO  m = 1, surf_usm_h%ns
+                       surf_usm_h%t_wall_av(iwl,m) =                           &
+                                          surf_usm_h%t_wall_av(iwl,m) /        &
+                                          REAL( average_count_3d, kind=wp )
+                    ENDDO
+                    DO  l = 0, 3
+                       DO  m = 1, surf_usm_v(l)%ns
+                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
+                                          surf_usm_v(l)%t_wall_av(iwl,m) /     &
+                                          REAL( average_count_3d, kind=wp )
+                       ENDDO
                     ENDDO
 
@@ -1493,9 +1795,8 @@
         TYPE(c_ptr)                                 :: lad_s_rma_p     !< allocated c pointer
         INTEGER(kind=MPI_ADDRESS_KIND)              :: size_lad_rma
-    
+!    
 !--     calculation of the SVF
         CALL location_message( '    calculation of SVF and CSF', .TRUE. )
-
-
+!
 !--     precalculate face areas for different face directions using normal vector
         DO d = 0, 9
@@ -1524,10 +1825,20 @@
 #if defined( __parallel )
         ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
-        nzterrl = nzb_s_inner(nys:nyn,nxl:nxr)
+        nzterrl = MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
+                               ), DIM = 1                                      &
+                        ) - 1  ! = nzb_s_inner(nys:nyn,nxl:nxr)
         CALL MPI_AllGather( nzterrl, nnx*nny, MPI_INTEGER, &
                             nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
         DEALLOCATE(nzterrl)
 #else
-        nzterr = RESHAPE( nzb_s_inner(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) )
+        nzterr = RESHAPE( MAXLOC(                                              &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
+                               ), DIM = 1                                      &
+                                ) - 1,                                         &
+                          (/(nx+1)*(ny+1)/)                                    &
+                        )
 #endif
         IF ( plant_canopy )  THEN
@@ -1579,5 +1890,10 @@
             DO i = nxl, nxr
                 DO j = nys, nyn
-                    k = nzb_s_inner(j, i)
+                    k = MAXLOC(                                                &
+                                MERGE( 1, 0,                                   &
+                                       BTEST( wall_flags_0(:,j,i), 12 )        &
+                                     ), DIM = 1                                &
+                              ) - 1
+
                     usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
                 ENDDO
@@ -1972,5 +2288,4 @@
         CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
        
-
     END SUBROUTINE usm_calc_svf
 
@@ -2089,5 +2404,6 @@
         INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
         INTEGER(iwp)                                           :: ids,isurf,isvf,isurfs,isurflt
-        INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat
+        INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
+        INTEGER(iwp)                                           ::  k_topo    !< topography top index
 
         dirstart = (/ startland, startwall, startwall, startwall, startwall /)
@@ -2139,42 +2455,71 @@
           CASE ( 'usm_surfz' )
 !--           array of lw radiation falling to local surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == ids )  THEN
-                      IF ( surfl(id,isurf) == iroof )  THEN
-                          temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) =             &
-                                  max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)),  &
-                                      REAL(surfl(iz,isurf),wp))
-                      ELSE
-                          temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) =             &
-                                  max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)),  &
-                                      REAL(surfl(iz,isurf),wp)+1.0_wp)
-                      ENDIF
-                  ENDIF
+              DO  m = 1, surf_usm_h%ns
+                 i = surf_usm_h%i(m)
+                 j = surf_usm_h%j(m)
+                 k = surf_usm_h%k(m)
+                 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
+              ENDDO
+              DO  l = 0, 3
+                 DO  m = 1, surf_usm_v(l)%ns
+                    i = surf_usm_v(l)%i(m)
+                    j = surf_usm_v(l)%j(m)
+                    k = surf_usm_v(l)%k(m)
+                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
+                 ENDDO
               ENDDO
 
           CASE ( 'usm_surfcat' )
 !--           surface category
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surface_types(isurf)
-                 ENDIF
+              DO  m = 1, surf_usm_h%ns
+                 i = surf_usm_h%i(m)
+                 j = surf_usm_h%j(m)
+                 k = surf_usm_h%k(m)
+                 temp_pf(k,j,i) = surf_usm_h%surface_types(m)
+              ENDDO
+              DO  l = 0, 3
+                 DO  m = 1, surf_usm_v(l)%ns
+                    i = surf_usm_v(l)%i(m)
+                    j = surf_usm_v(l)%j(m)
+                    k = surf_usm_v(l)%k(m)
+                    temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
+                 ENDDO
               ENDDO
               
           CASE ( 'usm_surfalb' )
 !--           surface albedo
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = albedo_surf(isurf)
-                 ENDIF
+              DO  m = 1, surf_usm_h%ns
+                 i = surf_usm_h%i(m)
+                 j = surf_usm_h%j(m)
+                 k = surf_usm_h%k(m)
+                 temp_pf(k,j,i) = surf_usm_h%albedo_surf(m)
+              ENDDO
+              DO  l = 0, 3
+                 DO  m = 1, surf_usm_v(l)%ns
+                    i = surf_usm_v(l)%i(m)
+                    j = surf_usm_v(l)%j(m)
+                    k = surf_usm_v(l)%k(m)
+                    temp_pf(k,j,i) = surf_usm_v(l)%albedo_surf(m)
+                 ENDDO
               ENDDO
               
           CASE ( 'usm_surfemis' )
 !--           surface albedo
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = emiss_surf(isurf)
-                 ENDIF
+              DO  m = 1, surf_usm_h%ns
+                 i = surf_usm_h%i(m)
+                 j = surf_usm_h%j(m)
+                 k = surf_usm_h%k(m)
+                 temp_pf(k,j,i) = surf_usm_h%emiss_surf(m)
               ENDDO
-              
+              DO  l = 0, 3
+                 DO  m = 1, surf_usm_v(l)%ns
+                    i = surf_usm_v(l)%i(m)
+                    j = surf_usm_v(l)%j(m)
+                    k = surf_usm_v(l)%k(m)
+                    temp_pf(k,j,i) = surf_usm_v(l)%emiss_surf(m)
+                 ENDDO
+              ENDDO
+!
+!-- Not adjusted so far              
           CASE ( 'usm_svf', 'usm_dif' )
 !--           shape view factors or iradiance factors to selected surface
@@ -2197,13 +2542,35 @@
           CASE ( 'usm_rad_net' )
 !--           array of complete radiation balance
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_l(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
+              IF ( av == 0 )  THEN
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
+                    ENDDO
+                 ENDDO
+              ELSE
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
+                    ENDDO
+                 ENDDO
+              ENDIF
 
           CASE ( 'usm_rad_insw' )
@@ -2268,16 +2635,4 @@
               ENDDO
 
-          CASE ( 'usm_rad_inlwdif' )
-!--           array of sw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
           CASE ( 'usm_rad_inlwref' )
 !--           array of lw radiation falling to surface from reflections
@@ -2339,65 +2694,175 @@
                  ENDIF
               ENDDO
-
+ 
           CASE ( 'usm_rad_hf' )
 !--           array of heat flux from radiation for surfaces after all reflections
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
+              IF ( av == 0 )  THEN
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%surfhf(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
+                    ENDDO
+                 ENDDO
+              ELSE
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
+                    ENDDO
+                 ENDDO
+              ENDIF
+ 
           CASE ( 'usm_wshf' )
 !--           array of sensible heat flux from surfaces
-!--           horizontal surfaces
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
+              IF ( av == 0 )  THEN
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
+                    ENDDO
+                 ENDDO
+              ELSE
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
+                    ENDDO
+                 ENDDO
+              ENDIF
+
 
           CASE ( 'usm_wghf' )
 !--           array of heat flux from ground (land, wall, roof)
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
+              IF ( av == 0 )  THEN
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
+                    ENDDO
+                 ENDDO
+              ELSE
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
+                    ENDDO
+                 ENDDO
+              ENDIF
 
           CASE ( 'usm_t_surf' )
 !--           surface temperature for surfaces
-              DO isurf = max(startenergy,dirstart(ids)), min(endenergy,dirend(ids))
-                 IF ( surfl(id,isurf) == ids )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
+              IF ( av == 0 )  THEN
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = t_surf_h(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = t_surf_v(l)%t(m)
+                    ENDDO
+                 ENDDO
+              ELSE
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m)
+                    ENDDO
+                 ENDDO
+              ENDIF
               
           CASE ( 'usm_t_wall' )
 !--           wall temperature for  iwl layer of walls and land
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall(iwl,isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall_av(iwl,isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
+              IF ( av == 0 )  THEN
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = t_wall_h(iwl,m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
+                    ENDDO
+                 ENDDO
+              ELSE
+                 DO  m = 1, surf_usm_h%ns
+                    i = surf_usm_h%i(m)
+                    j = surf_usm_h%j(m)
+                    k = surf_usm_h%k(m)
+                    temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
+                 ENDDO
+                 DO  l = 0, 3
+                    DO  m = 1, surf_usm_v(l)%ns
+                       i = surf_usm_v(l)%i(m)
+                       j = surf_usm_v(l)%j(m)
+                       k = surf_usm_v(l)%k(m)
+                       temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
+                    ENDDO
+                 ENDDO
+              ENDIF
              
           CASE DEFAULT
@@ -2408,4 +2873,6 @@
 !--     fill out array local_pf which is subsequently treated by data_output_3d
         CALL exchange_horiz( temp_pf, nbgp )
+!
+!--  To Do: why reversed loop order
         DO j = nysg,nyng
             DO i = nxlg,nxrg
@@ -2589,5 +3056,5 @@
         IMPLICIT NONE
 
-        INTEGER(iwp) ::  k, l            !< running indices
+        INTEGER(iwp) ::  k, l, m            !< running indices
         
         CALL location_message( '    initialization of wall surface model', .TRUE. )
@@ -2596,26 +3063,55 @@
 !--     Temperature is defined at the center of the wall layers,
 !--     whereas gradients/fluxes are defined at the edges (_stag)
-        DO l = nzb_wall, nzt_wall
-           zwn(l) = zwn_default(l)
+        DO k = nzb_wall, nzt_wall
+           zwn(k) = zwn_default(k)
         ENDDO
-        
-!--     apply for all particular wall grids
-        DO l = startenergy, endenergy
-           zw(:,l) = zwn(:) * thickness_wall(l)
-           dz_wall(nzb_wall,l) = zw(nzb_wall,l)
+!        
+!--     apply for all particular surface grids. First for horizontal surfaces
+        DO  m = 1, surf_usm_h%ns
+           surf_usm_h%zw(:,m)             = zwn(:) *                           &
+                                            surf_usm_h%thickness_wall(m)
+           surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
            DO k = nzb_wall+1, nzt_wall
-               dz_wall(k,l) = zw(k,l) - zw(k-1,l)
+               surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) -                  &
+                                         surf_usm_h%zw(k-1,m)
            ENDDO
            
-           dz_wall(nzt_wall+1,l) = dz_wall(nzt_wall,l)
+           surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
 
            DO k = nzb_wall, nzt_wall-1
-               dz_wall_stag(k,l) = 0.5 * (dz_wall(k+1,l) + dz_wall(k,l))
+               surf_usm_h%dz_wall_stag(k,m) = 0.5 * (                          &
+                           surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
            ENDDO
-           dz_wall_stag(nzt_wall,l) = dz_wall(nzt_wall,l)
+           surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
         ENDDO
-        
-        ddz_wall      = 1.0_wp / dz_wall
-        ddz_wall_stag = 1.0_wp / dz_wall_stag
+        surf_usm_h%ddz_wall      = 1.0_wp / surf_usm_h%dz_wall
+        surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag
+!        
+!--     For vertical surfaces
+        DO  l = 0, 3
+           DO  m = 1, surf_usm_v(l)%ns
+              surf_usm_v(l)%zw(:,m)             = zwn(:) *                     &
+                                                  surf_usm_v(l)%thickness_wall(m)
+              surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
+              DO k = nzb_wall+1, nzt_wall
+                  surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) -         &
+                                               surf_usm_v(l)%zw(k-1,m)
+              ENDDO
+           
+              surf_usm_v(l)%dz_wall(nzt_wall+1,m) =                            &
+                                              surf_usm_v(l)%dz_wall(nzt_wall,m)
+
+              DO k = nzb_wall, nzt_wall-1
+                  surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * (                    &
+                                                surf_usm_v(l)%dz_wall(k+1,m) + &
+                                                surf_usm_v(l)%dz_wall(k,m) )
+              ENDDO
+              surf_usm_v(l)%dz_wall_stag(nzt_wall,m) =                         &
+                                              surf_usm_v(l)%dz_wall(nzt_wall,m)
+           ENDDO
+           surf_usm_v(l)%ddz_wall      = 1.0_wp / surf_usm_v(l)%dz_wall
+           surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag
+        ENDDO      
+
         
         CALL location_message( '    wall structures filed out', .TRUE. )
@@ -2635,5 +3131,5 @@
         IMPLICIT NONE
 
-        INTEGER(iwp) ::  i, j, k, l            !< running indices
+        INTEGER(iwp) ::  i, j, k, l, m            !< running indices
         REAL(wp)     ::  c, d, tin, exn
         
@@ -2699,42 +3195,83 @@
         ELSE
        
-!--         Calculate initial surface temperature
+!--         Calculate initial surface temperature from pt of adjacent gridbox
             exn = ( surface_pressure / 1000.0_wp )**0.286_wp
 
-            DO l = startenergy, endenergy
-                k = surfl(iz,l)
-                j = surfl(iy,l)
-                i = surfl(ix,l)
-
-!--              Initial surface temperature set from pt of adjacent gridbox
-                t_surf(l) = pt(k,j,i) * exn
+!
+!--         At horizontal surfaces. Please note, t_surf_h is defined on a 
+!--         different data type, but with the same dimension.
+            DO  m = 1, surf_usm_h%ns
+               i = surf_usm_h%i(m)            
+               j = surf_usm_h%j(m)
+               k = surf_usm_h%k(m)
+
+               t_surf_h(m) = pt(k,j,i) * exn
             ENDDO
+!
+!--         At vertical surfaces.
+            DO  l = 0, 3
+               DO  m = 1, surf_usm_v(l)%ns
+                  i = surf_usm_v(l)%i(m)            
+                  j = surf_usm_v(l)%j(m)
+                  k = surf_usm_v(l)%k(m)
+
+                  t_surf_v(l)%t(m) = pt(k,j,i) * exn
+               ENDDO
+            ENDDO
+
       
 !--         initial values for t_wall
 !--         outer value is set to surface temperature
 !--         inner value is set to wall_inner_temperature
-!--         and profile is logaritmic (linear in nz)
-            DO l = startenergy, endenergy
-                IF ( isroof_surf(l) )  THEN
-                    tin = roof_inner_temperature
-                ELSE IF ( surf(id,l)==iroof )  THEN
-                    tin = soil_inner_temperature
-                ELSE
-                    tin = wall_inner_temperature
-                ENDIF
-                DO k = nzb_wall, nzt_wall+1
-                    c = REAL(k-nzb_wall,wp)/REAL(nzt_wall+1-nzb_wall,wp)
-                    t_wall(k,:) = (1.0_wp-c)*t_surf(:) + c*tin
-                ENDDO
+!--         and profile is logaritmic (linear in nz).
+!--         Horizontal surfaces
+            DO  m = 1, surf_usm_h%ns
+!
+!--            Roof
+               IF ( surf_usm_h%isroof_surf(m) )  THEN
+                   tin = roof_inner_temperature
+!
+!--            Normal land surface
+               ELSE 
+                   tin = soil_inner_temperature
+               ENDIF
+
+               DO k = nzb_wall, nzt_wall+1
+                   c = REAL( k - nzb_wall, wp ) /                              &
+                       REAL( nzt_wall + 1 - nzb_wall , wp )
+
+                   t_wall_h(k,m) = ( 1.0_wp - c ) * t_surf_h(m) + c * tin
+               ENDDO
             ENDDO
+!
+!--         Vertical surfaces
+            DO  l = 0, 3
+               DO  m = 1, surf_usm_v(l)%ns
+!
+!--               Inner wall
+                  tin = wall_inner_temperature
+
+                  DO k = nzb_wall, nzt_wall+1
+                     c = REAL( k - nzb_wall, wp ) /                            &
+                         REAL( nzt_wall + 1 - nzb_wall , wp )
+
+                     t_wall_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_v(l)%t(m) +  &
+                                          c * tin
+                  ENDDO
+               ENDDO
+            ENDDO
+
         ENDIF
         
 !--    
-!--        Possibly DO user-defined actions (e.g. define heterogeneous wall surface)
+!--     Possibly DO user-defined actions (e.g. define heterogeneous wall surface)
         CALL user_init_urban_surface
 
 !--     initialize prognostic values for the first timestep
-        t_surf_p = t_surf
-        t_wall_p = t_wall
+        t_surf_h_p = t_surf_h
+        t_surf_v_p = t_surf_v
+
+        t_wall_h_p = t_wall_h
+        t_wall_v_p = t_wall_v
         
 !--     Adjust radiative fluxes for urban surface at model start
@@ -2759,54 +3296,113 @@
         IMPLICIT NONE
 
-        INTEGER(iwp) ::  i,j,k,l,kw                      !< running indices
+        INTEGER(iwp) ::  i,j,k,l,kw, m                      !< running indices
 
         REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend  !< tendency
 
-                                               
-        DO l = startenergy, endenergy
-!--         calculate frequently used parameters
-            k = surfl(iz,l)
-            j = surfl(iy,l)
-            i = surfl(ix,l)
-
-            !
-!--         prognostic equation for ground/wall/roof temperature t_wall
-            wtend(:) = 0.0_wp
-            wtend(nzb_wall) = (1.0_wp/rho_c_wall(nzb_wall,l)) *                     &
-                       ( lambda_h(nzb_wall,l) * ( t_wall(nzb_wall+1,l)              &
-                         - t_wall(nzb_wall,l) ) * ddz_wall(nzb_wall+1,l)            &
-                         + wghf_eb(l) ) * ddz_wall_stag(nzb_wall,l)
+!
+!--     For horizontal surfaces                                   
+        DO  m = 1, surf_usm_h%ns
+!
+!--        Obtain indices
+           i = surf_usm_h%i(m)            
+           j = surf_usm_h%j(m)
+           k = surf_usm_h%k(m)
+!
+!--        prognostic equation for ground/roof temperature t_wall_h
+           wtend(:) = 0.0_wp
+           wtend(nzb_wall) = (1.0_wp / surf_usm_h%rho_c_wall(nzb_wall,m)) *    &
+                                      ( surf_usm_h%lambda_h(nzb_wall,m) *      &
+                                        ( t_wall_h(nzb_wall+1,m)               &
+                                        - t_wall_h(nzb_wall,m) ) *             &
+                                        surf_usm_h%ddz_wall(nzb_wall+1,m)      &
+                                      + surf_usm_h%wghf_eb(m) ) *              &
+                                        surf_usm_h%ddz_wall_stag(nzb_wall,m)
             
-            DO  kw = nzb_wall+1, nzt_wall
-                wtend(kw) = (1.0_wp/rho_c_wall(kw,l))                               &
-                              * (   lambda_h(kw,l)                                  &
-                                 * ( t_wall(kw+1,l) - t_wall(kw,l) )                &
-                                 * ddz_wall(kw+1,l)                                 &
-                              - lambda_h(kw-1,l)                                    &
-                                 * ( t_wall(kw,l) - t_wall(kw-1,l) )                &
-                                 * ddz_wall(kw,l)                                   &
-                              ) * ddz_wall_stag(kw,l)
+           DO  kw = nzb_wall+1, nzt_wall
+               wtend(kw) = (1.0_wp / surf_usm_h%rho_c_wall(kw,m))              &
+                              * (   surf_usm_h%lambda_h(kw,m)                  &
+                                 * ( t_wall_h(kw+1,m) - t_wall_h(kw,m) )       &
+                                 * surf_usm_h%ddz_wall(kw+1,m)                 &
+                              - surf_usm_h%lambda_h(kw-1,m)                    &
+                                 * ( t_wall_h(kw,m) - t_wall_h(kw-1,m) )       &
+                                 * surf_usm_h%ddz_wall(kw,m)                   &
+                              ) * surf_usm_h%ddz_wall_stag(kw,m)
             ENDDO
 
-            t_wall_p(nzb_wall:nzt_wall,l) = t_wall(nzb_wall:nzt_wall,l)             &
-                                             + dt_3d * ( tsc(2)                     &
-                                             * wtend(nzb_wall:nzt_wall) + tsc(3)    &
-                                             * tt_wall_m(nzb_wall:nzt_wall,l) )   
+           t_wall_h_p(nzb_wall:nzt_wall,m) = t_wall_h(nzb_wall:nzt_wall,m)     &
+                                 + dt_3d * ( tsc(2)                            &
+                                 * wtend(nzb_wall:nzt_wall) + tsc(3)           &
+                                 * surf_usm_h%tt_wall_m(nzb_wall:nzt_wall,m) )   
             
-            !
-!--         calculate t_wall tendencies for the next Runge-Kutta step
-            IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
+!
+!--        calculate t_wall tendencies for the next Runge-Kutta step
+           IF ( timestep_scheme(1:5) == 'runge' )  THEN
+               IF ( intermediate_timestep_count == 1 )  THEN
+                  DO  kw = nzb_wall, nzt_wall
+                     surf_usm_h%tt_wall_m(kw,m) = wtend(kw)
+                  ENDDO
+               ELSEIF ( intermediate_timestep_count <                          &
+                        intermediate_timestep_count_max )  THEN
                    DO  kw = nzb_wall, nzt_wall
-                      tt_wall_m(kw,l) = wtend(kw)
+                      surf_usm_h%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) +    &
+                                         5.3125_wp * surf_usm_h%tt_wall_m(kw,m)
                    ENDDO
-                ELSEIF ( intermediate_timestep_count <                              &
-                         intermediate_timestep_count_max )  THEN
-                    DO  kw = nzb_wall, nzt_wall
-                        tt_wall_m(kw,l) = -9.5625_wp * wtend(kw) + 5.3125_wp        &
-                                         * tt_wall_m(kw,l)
-                    ENDDO
-                ENDIF
-            ENDIF
+               ENDIF
+           ENDIF
+        ENDDO
+!
+!--     For vertical surfaces     
+        DO  l = 0, 3                              
+           DO  m = 1, surf_usm_v(l)%ns
+!
+!--           Obtain indices
+              i = surf_usm_v(l)%i(m)            
+              j = surf_usm_v(l)%j(m)
+              k = surf_usm_v(l)%k(m)
+!
+!--           prognostic equation for wall temperature t_wall_v
+              wtend(:) = 0.0_wp
+              wtend(nzb_wall) = (1.0_wp / surf_usm_v(l)%rho_c_wall(nzb_wall,m)) * &
+                                      ( surf_usm_v(l)%lambda_h(nzb_wall,m) *   &
+                                        ( t_wall_v(l)%t(nzb_wall+1,m)          &
+                                        - t_wall_v(l)%t(nzb_wall,m) ) *        &
+                                        surf_usm_v(l)%ddz_wall(nzb_wall+1,m)   &
+                                      + surf_usm_v(l)%wghf_eb(m) ) *           &
+                                        surf_usm_v(l)%ddz_wall_stag(nzb_wall,m)
+            
+              DO  kw = nzb_wall+1, nzt_wall
+                  wtend(kw) = (1.0_wp / surf_usm_v(l)%rho_c_wall(kw,m))        &
+                           * (   surf_usm_v(l)%lambda_h(kw,m)                  &
+                              * ( t_wall_v(l)%t(kw+1,m) - t_wall_v(l)%t(kw,m) )&
+                              * surf_usm_v(l)%ddz_wall(kw+1,m)                 &
+                           - surf_usm_v(l)%lambda_h(kw-1,m)                    &
+                              * ( t_wall_v(l)%t(kw,m) - t_wall_v(l)%t(kw-1,m) )&
+                              * surf_usm_v(l)%ddz_wall(kw,m)                   &
+                              ) * surf_usm_v(l)%ddz_wall_stag(kw,m)
+               ENDDO
+
+              t_wall_v_p(l)%t(nzb_wall:nzt_wall,m) =                           &
+                                   t_wall_v(l)%t(nzb_wall:nzt_wall,m)          &
+                                 + dt_3d * ( tsc(2)                            &
+                                 * wtend(nzb_wall:nzt_wall) + tsc(3)           &
+                                 * surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall,m) )   
+            
+!
+!--           calculate t_wall tendencies for the next Runge-Kutta step
+              IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                  IF ( intermediate_timestep_count == 1 )  THEN
+                     DO  kw = nzb_wall, nzt_wall
+                        surf_usm_v(l)%tt_wall_m(kw,m) = wtend(kw)
+                     ENDDO
+                  ELSEIF ( intermediate_timestep_count <                       &
+                           intermediate_timestep_count_max )  THEN
+                      DO  kw = nzb_wall, nzt_wall
+                         surf_usm_v(l)%tt_wall_m(kw,m) =                       &
+                                     - 9.5625_wp * wtend(kw) +                 &
+                                       5.3125_wp * surf_usm_v(l)%tt_wall_m(kw,m)
+                      ENDDO
+                  ENDIF
+              ENDIF
+           ENDDO
         ENDDO
 
@@ -2857,9 +3453,7 @@
 !--    Read user-defined namelist
        READ ( 11, urban_surface_par )
-
 !
 !--    Set flag that indicates that the land surface model is switched on
        urban_surface = .TRUE.
-       
 
  10    CONTINUE
@@ -2879,5 +3473,5 @@
         IMPLICIT NONE
        
-        INTEGER(iwp)               :: i, j, k, kk, is, js, d, ku, refstep
+        INTEGER(iwp)               :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
         INTEGER(iwp)               :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb
         INTEGER(iwp), DIMENSION(4) :: bdycross
@@ -2892,4 +3486,5 @@
         REAL(wp)                   :: pc_box_area, pc_abs_frac, pc_abs_eff
         INTEGER(iwp)               :: pc_box_dimshift !< transform for best accuracy
+        INTEGER(iwp), DIMENSION(0:3) :: reorder = (/ 1, 0, 3, 2 /)
         
         
@@ -2949,23 +3544,55 @@
 !--     First pass: direct + diffuse irradiance
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        surfinswdir   = 0._wp
-        surfinswdif   = 0._wp
-        surfinlwdif   = 0._wp
-        surfins   = 0._wp
-        surfinl   = 0._wp
-        surfoutsl    = 0._wp
-        surfoutll    = 0._wp
+        surfinswdir   = 0._wp !nsurfl
+        surfinswdif   = 0._wp !nsurfl
+        surfinlwdif   = 0._wp !nsurfl
+        surfins       = 0._wp !nsurfl
+        surfinl       = 0._wp !nsurfl
+        surfoutsl(:)  = 0.0_wp !start-end
+        surfoutll(:)  = 0.0_wp !start-end
         
 !--     Set up thermal radiation from surfaces
 !--     emiss_surf is defined only for surfaces for which energy balance is calculated
-        surfoutll(startenergy:endenergy) = emiss_surf(startenergy:endenergy) * sigma_sb   &
-                                           * t_surf(startenergy:endenergy)**4
+!--     Workaround: reorder surface data type back on 1D array including all surfaces, 
+!--     which implies to reorder horizontal and vertical surfaces
+!
+!--     Horizontal walls
+        mm = 1
+        DO  i = nxl, nxr
+           DO  j = nys, nyn
+
+              DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+                 surfoutll(mm) = surf_usm_h%emiss_surf(m) * sigma_sb   &
+                                     * t_surf_h(m)**4
+                 albedo_surf(mm) = surf_usm_h%albedo_surf(m)
+                 emiss_surf(mm)  = surf_usm_h%emiss_surf(m)
+                 mm = mm + 1
+              ENDDO
+           ENDDO
+        ENDDO
+!
+!--     Vertical walls
+        DO  i = nxl, nxr
+           DO  j = nys, nyn
+              DO  ll = 0, 3
+                 l = reorder(ll)
+                 DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
+                    surfoutll(mm) = surf_usm_v(l)%emiss_surf(m) * sigma_sb   &
+                                     * t_surf_v(l)%t(m)**4
+                    albedo_surf(mm) = surf_usm_v(l)%albedo_surf(m)
+                    emiss_surf(mm) = surf_usm_v(l)%emiss_surf(m)
+                    mm = mm + 1
+                 ENDDO
+              ENDDO
+           ENDDO
+        ENDDO
         
 #if defined( __parallel )
 !--     might be optimized and gather only values relevant for current processor
+        
         CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, &
-                            surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
+                            surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
 #else
-        surfoutl(:) = surfoutll(:)
+        surfoutl(:) = surfoutll(:) !nsurf global
 #endif
         
@@ -3005,8 +3632,8 @@
                 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
             ENDIF
-            
-            IF ( zenith(0) > 0  .AND.  all( surf(:, isurfsrc) == bdycross ) )  THEN
+
+            IF ( zenith(0) > 0  .AND.  all( surf(1:4,isurfsrc) == bdycross ) )  THEN
 !--             found svf between model boundary and the face => face isn't shaded
-                surfinswdir(isurf) = rad_sw_in_dir(j, i) &
+                surfinswdir(isurf) = rad_sw_in_dir(j,i) &
                     * costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0)
 
@@ -3050,5 +3677,5 @@
                 ENDIF
 
-                IF ( zenith(0) > 0  .AND.  all( surf(:, isurfsrc) == bdycross ) )  THEN
+                IF ( zenith(0) > 0  .AND.  all( surf(1:4,isurfsrc) == bdycross ) )  THEN
 !--                 found svf between model boundary and the pcgb => pcgb isn't shaded
                     pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
@@ -3058,4 +3685,5 @@
             ENDDO
         ENDIF
+
         surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy)
         surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy)
@@ -3064,6 +3692,6 @@
         surfoutsw(:) = 0.0_wp
         surfoutlw(:) = surfoutll(:)
-        surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
-                                      - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
+!         surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
+!                                       - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
         
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -3117,6 +3745,6 @@
             surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy)
             surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy)
-            surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
-                                          - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
+!             surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
+!                                           - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
         
         ENDDO
@@ -3129,9 +3757,72 @@
                 i = pcbl(ix, ipcgb)
                 k = pcbl(iz, ipcgb)
-                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
+                kk = k - MAXLOC(                                               &
+                                MERGE( 1, 0,                                   &
+                                       BTEST( wall_flags_0(:,j,i), 12 )        &
+                                     ), DIM = 1                                &
+                               ) - 1 ! kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
                 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
                     * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
             ENDDO
         ENDIF
+!
+!--     Transfer radiation arrays required for energy balance to the respective data types
+        DO  i = startenergy, endenergy
+           m  = surfl(5,i)          
+!
+!--        upward-facing
+           IF ( surfl(1,i) == 0 )  THEN
+              surf_usm_h%rad_in_sw(m)  = surfinsw(i) 
+              surf_usm_h%rad_out_sw(m) = surfoutsw(i) 
+              surf_usm_h%rad_in_lw(m)  = surfinlw(i)
+              surf_usm_h%rad_out_lw(m) = surfoutlw(i)
+!
+!--        southward-facding
+           ELSEIF ( surfl(1,i) == 1 )  THEN
+              surf_usm_v(1)%rad_in_sw(m)  = surfinsw(i) 
+              surf_usm_v(1)%rad_out_sw(m) = surfoutsw(i) 
+              surf_usm_v(1)%rad_in_lw(m)  = surfinlw(i)
+              surf_usm_v(1)%rad_out_lw(m) = surfoutlw(i)
+!
+!--        northward-facding
+           ELSEIF ( surfl(1,i) == 2 )  THEN
+              surf_usm_v(0)%rad_in_sw(m)  = surfinsw(i) 
+              surf_usm_v(0)%rad_out_sw(m) = surfoutsw(i) 
+              surf_usm_v(0)%rad_in_lw(m)  = surfinlw(i)
+              surf_usm_v(0)%rad_out_lw(m) = surfoutlw(i)
+!
+!--        westward-facding
+           ELSEIF ( surfl(1,i) == 3 )  THEN
+              surf_usm_v(3)%rad_in_sw(m)  = surfinsw(i) 
+              surf_usm_v(3)%rad_out_sw(m) = surfoutsw(i) 
+              surf_usm_v(3)%rad_in_lw(m)  = surfinlw(i)
+              surf_usm_v(3)%rad_out_lw(m) = surfoutlw(i)
+!
+!--        eastward-facing
+           ELSEIF ( surfl(1,i) == 4 )  THEN
+              surf_usm_v(2)%rad_in_sw(m)  = surfinsw(i) 
+              surf_usm_v(2)%rad_out_sw(m) = surfoutsw(i) 
+              surf_usm_v(2)%rad_in_lw(m)  = surfinlw(i)
+              surf_usm_v(2)%rad_out_lw(m) = surfoutlw(i)
+           ENDIF
+
+        ENDDO
+
+
+        DO  m = 1, surf_usm_h%ns
+           surf_usm_h%surfhf(m) = surf_usm_h%rad_in_sw(m)  +                   &
+                                  surf_usm_h%rad_in_lw(m)  -                   &
+                                  surf_usm_h%rad_out_sw(m) -                   &
+                                  surf_usm_h%rad_out_lw(m)
+        ENDDO
+
+        DO  l = 0, 3
+           DO  m = 1, surf_usm_v(l)%ns
+              surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_in_sw(m)  +          &
+                                        surf_usm_v(l)%rad_in_lw(m)  -          &
+                                        surf_usm_v(l)%rad_out_sw(m) -          &
+                                        surf_usm_v(l)%rad_out_lw(m)
+           ENDDO
+        ENDDO
 
 !--     return surface radiation to horizontal surfaces 
@@ -3207,5 +3898,5 @@
         REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
 
-
+!
 !--     Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
 !--     the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
@@ -3351,5 +4042,4 @@
         visible = .TRUE.
 
-        
     END SUBROUTINE usm_raytrace
     
@@ -3470,24 +4160,103 @@
                 SELECT CASE ( TRIM( variable_chr ) )
                 
-                   CASE ( 't_surf' )
+                   CASE ( 't_surf_h' )
 #if defined( __nopointer )                   
-                      IF ( .NOT.  ALLOCATED( t_surf ) )                         &
-                         ALLOCATE( t_surf(startenergy:endenergy) )
-                      READ ( 13 )  t_surf
+                      IF ( .NOT.  ALLOCATED( t_surf_h ) )                      &
+                         ALLOCATE( t_surf_h(1:surf_usm_h%ns) )
+                      READ ( 13 )  t_surf_h
 #else                      
-                      IF ( .NOT.  ALLOCATED( t_surf_1 ) )                         &
-                         ALLOCATE( t_surf_1(startenergy:endenergy) )
-                      READ ( 13 )  t_surf_1
+                      IF ( .NOT.  ALLOCATED( t_surf_h_1 ) )                    &
+                         ALLOCATE( t_surf_h_1(1:surf_usm_h%ns) )
+                      READ ( 13 )  t_surf_h_1
 #endif
-
-                   CASE ( 't_wall' )
+                   CASE ( 't_surf_v(0)' )
+#if defined( __nopointer )                   
+                      IF ( .NOT.  ALLOCATED( t_surf_v(0)%t ) )                 &
+                         ALLOCATE( t_surf_v(0)%t(1:surf_usm_v(0)%ns) )
+                      READ ( 13 )  t_surf_v(0)%t
+#else                      
+                      IF ( .NOT.  ALLOCATED( t_surf_v_1(0)%t ) )               &
+                         ALLOCATE( t_surf_v_1(0)%t(1:surf_usm_v(0)%ns) )
+                      READ ( 13 )  t_surf_v_1(0)%t
+#endif
+                   CASE ( 't_surf_v(1)' )
+#if defined( __nopointer )                   
+                      IF ( .NOT.  ALLOCATED( t_surf_v(1)%t ) )                 &
+                         ALLOCATE( t_surf_v(1)%t(1:surf_usm_v(1)%ns) )
+                      READ ( 13 )  t_surf_v(1)%t
+#else                      
+                      IF ( .NOT.  ALLOCATED( t_surf_v_1(1)%t ) )               &
+                         ALLOCATE( t_surf_v_1(1)%t(1:surf_usm_v(1)%ns) )
+                      READ ( 13 )  t_surf_v_1(1)%t
+#endif
+                   CASE ( 't_surf_v(2)' )
+#if defined( __nopointer )                   
+                      IF ( .NOT.  ALLOCATED( t_surf_v(2)%t ) )                 &
+                         ALLOCATE( t_surf_v(2)%t(1:surf_usm_v(2)%ns) )
+                      READ ( 13 )  t_surf_v(2)%t
+#else                      
+                      IF ( .NOT.  ALLOCATED( t_surf_v_1(2)%t ) )               &
+                         ALLOCATE( t_surf_v_1(2)%t(1:surf_usm_v(2)%ns) )
+                      READ ( 13 )  t_surf_v_1(2)%t
+#endif
+                   CASE ( 't_surf_v(3)' )
+#if defined( __nopointer )                   
+                      IF ( .NOT.  ALLOCATED( t_surf_v(3)%t ) )                 &
+                         ALLOCATE( t_surf_v(3)%t(1:surf_usm_v(3)%ns) )
+                      READ ( 13 )  t_surf_v(3)%t
+#else                      
+                      IF ( .NOT.  ALLOCATED( t_surf_v_1(3)%t ) )               &
+                         ALLOCATE( t_surf_v_1(3)%t(1:surf_usm_v(3)%ns) )
+                      READ ( 13 )  t_surf_v_1(3)%t
+#endif
+                   CASE ( 't_wall_h' )
 #if defined( __nopointer )
-                      IF ( .NOT.  ALLOCATED( t_wall ) )                         &
-                         ALLOCATE( t_wall(nzb_wall:nzt_wall+1,startenergy:endenergy) )
-                      READ ( 13 )  t_wall
+                      IF ( .NOT.  ALLOCATED( t_wall_h ) )                      &
+                         ALLOCATE( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
+                      READ ( 13 )  t_wall_h
 #else
-                      IF ( .NOT.  ALLOCATED( t_wall_1 ) )                         &
-                         ALLOCATE( t_wall_1(nzb_wall:nzt_wall+1,startenergy:endenergy) )
-                      READ ( 13 )  t_wall_1
+                      IF ( .NOT.  ALLOCATED( t_wall_h_1 ) )                    &
+                         ALLOCATE( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
+                      READ ( 13 )  t_wall_h_1
+#endif
+                   CASE ( 't_wall_v(0)' )
+#if defined( __nopointer )
+                      IF ( .NOT.  ALLOCATED( t_wall_v(0)%t ) )                      &
+                         ALLOCATE( t_wall_v(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
+                      READ ( 13 )  t_wall_v(0)%t
+#else
+                      IF ( .NOT.  ALLOCATED( t_wall_v_1(0)%t ) )                    &
+                         ALLOCATE( t_wall_v_1(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
+                      READ ( 13 )  t_wall_v_1(0)%t
+#endif
+                   CASE ( 't_wall_v(1)' )
+#if defined( __nopointer )
+                      IF ( .NOT.  ALLOCATED( t_wall_v(1)%t ) )                      &
+                         ALLOCATE( t_wall_v(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
+                      READ ( 13 )  t_wall_v(1)%t
+#else
+                      IF ( .NOT.  ALLOCATED( t_wall_v_1(0)%t ) )                    &
+                         ALLOCATE( t_wall_v_1(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
+                      READ ( 13 )  t_wall_v_1(1)%t
+#endif
+                   CASE ( 't_wall_v(2)' )
+#if defined( __nopointer )
+                      IF ( .NOT.  ALLOCATED( t_wall_v(2)%t ) )                      &
+                         ALLOCATE( t_wall_v(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
+                      READ ( 13 )  t_wall_v(2)%t
+#else
+                      IF ( .NOT.  ALLOCATED( t_wall_v_1(2)%t ) )                    &
+                         ALLOCATE( t_wall_v_1(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
+                      READ ( 13 )  t_wall_v_1(2)%t
+#endif
+                   CASE ( 't_wall_v(3)' )
+#if defined( __nopointer )
+                      IF ( .NOT.  ALLOCATED( t_wall_v(3)%t ) )                      &
+                         ALLOCATE( t_wall_v(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
+                      READ ( 13 )  t_wall_v(3)%t
+#else
+                      IF ( .NOT.  ALLOCATED( t_wall_v_1(3)%t ) )                    &
+                         ALLOCATE( t_wall_v_1(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
+                      READ ( 13 )  t_wall_v_1(3)%t
 #endif
 
@@ -3593,5 +4362,5 @@
         INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg)   :: usm_par
         REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg)       :: usm_val
-        INTEGER(iwp)                                          :: k, l, d, iw, jw, kw, it, ip, ii, ij
+        INTEGER(iwp)                                          :: k, l, d, iw, jw, kw, it, ip, ii, ij, m
         INTEGER(iwp)                                          :: i, j
         INTEGER(iwp)                                          :: nz, roof, dirwe, dirsn
@@ -3733,152 +4502,198 @@
             ENDDO
         ENDDO
-        
-!--     assign the surface types to local surface array
-        DO  l = startenergy, endenergy
-            
-            d = surfl(id,l)
-            kw = surfl(iz,l)
-            j = surfl(iy,l)
-            i = surfl(ix,l)
-            IF ( d == iroof )  THEN
-!--             horizontal surface - land or roof
-                iw = i
-                jw = j
-                IF ( usm_par(5,jw,iw) == 0 )  THEN
-                    IF ( zu(kw) >= roof_height_limit )  THEN
-                        isroof_surf(l) = .TRUE.
-                        surface_types(l) = roof_category         !< default category for root surface
-                    ELSE
-                        isroof_surf(l) = .FALSE.
-                        surface_types(l) = land_category         !< default category for land surface
-                    ENDIF
-                    albedo_surf(l) = -1.0_wp
-                    thickness_wall(l) = -1.0_wp
-                ELSE
-                    IF ( usm_par(2,jw,iw)==0 )  THEN
-                        isroof_surf(l) = .FALSE.
-                        thickness_wall(l) = -1.0_wp
-                    ELSE
-                        isroof_surf(l) = .TRUE.
-                        thickness_wall(l) = usm_val(2,jw,iw)
-                    ENDIF
-                    surface_types(l) = usm_par(5,jw,iw)
-                    albedo_surf(l) = usm_val(1,jw,iw)
-                ENDIF
-            ELSE
-                SELECT CASE (d)
-                    CASE (iwest)
-                        iw = i
-                        jw = j
-                        ii = 6
-                        ij = 3
-                    CASE (ieast)
-                        iw = i-1
-                        jw = j
-                        ii = 6
-                        ij = 3
-                    CASE (isouth)
-                        iw = i
-                        jw = j
-                        ii = 12
-                        ij = 9
-                    CASE (inorth)
-                        iw = i
-                        jw = j-1
-                        ii = 12
-                        ij = 9
-                END SELECT
-                
-                IF ( kw <= usm_par(ii,jw,iw) )  THEN
-!--                 pedestrant zone
-                    isroof_surf(l) = .FALSE.
-                    IF ( usm_par(ii+1,jw,iw) == 0 )  THEN
-                        surface_types(l) = pedestrant_category   !< default category for wall surface in pedestrant zone
-                        albedo_surf(l) = -1.0_wp
-                        thickness_wall(l) = -1.0_wp
-                    ELSE
-                        surface_types(l) = usm_par(ii+1,jw,iw)
-                        albedo_surf(l) = usm_val(ij,jw,iw)
-                        thickness_wall(l) = usm_val(ij+1,jw,iw)
-                    ENDIF
-                ELSE IF ( kw <= usm_par(ii+2,jw,iw) )  THEN
-!--                 wall zone
-                    isroof_surf(l) = .FALSE.
-                    IF ( usm_par(ii+3,jw,iw) == 0 )  THEN
-                        surface_types(l) = wall_category         !< default category for wall surface
-                        albedo_surf(l) = -1.0_wp
-                        thickness_wall(l) = -1.0_wp
-                    ELSE
-                        surface_types(l) = usm_par(ii+3,jw,iw)
-                        albedo_surf(l) = usm_val(ij+2,jw,iw)
-                        thickness_wall(l) = usm_val(ij+3,jw,iw)
-                    ENDIF
-                ELSE IF ( kw <= usm_par(ii+4,jw,iw) )  THEN
-!--                 roof zone
-                    isroof_surf(l) = .TRUE.
-                    IF ( usm_par(ii+5,jw,iw) == 0 )  THEN
-                        surface_types(l) = roof_category         !< default category for roof surface
-                        albedo_surf(l) = -1.0_wp
-                        thickness_wall(l) = -1.0_wp
-                    ELSE
-                        surface_types(l) = usm_par(ii+5,jw,iw)
-                        albedo_surf(l) = usm_val(ij+4,jw,iw)
-                        thickness_wall(l) = usm_val(ij+5,jw,iw)
-                    ENDIF
-                ELSE
-!--                 something wrong
-                    CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 )
-                ENDIF
-            ENDIF
-            
-!--         find the type position
-            it = surface_types(l)
-            ip = -99999
-            DO k = 1, n_surface_types
-                IF ( surface_type_codes(k) == it )  THEN
+!        
+!--     Assign the surface types to the respective data type. 
+!--     First, for horizontal upward-facing surfaces. 
+        DO  m = 1, surf_usm_h%ns
+           iw = surf_usm_h%i(m)
+           jw = surf_usm_h%j(m)
+           kw = surf_usm_h%k(m)
+
+           IF ( usm_par(5,jw,iw) == 0 )  THEN
+              IF ( zu(kw) >= roof_height_limit )  THEN
+                 surf_usm_h%isroof_surf(m)   = .TRUE.
+                 surf_usm_h%surface_types(m) = roof_category         !< default category for root surface
+              ELSE
+                 surf_usm_h%isroof_surf(m)   = .FALSE.
+                 surf_usm_h%surface_types(m) = land_category         !< default category for land surface
+              ENDIF
+              surf_usm_h%albedo_surf(m)    = -1.0_wp
+              surf_usm_h%thickness_wall(m) = -1.0_wp
+           ELSE
+              IF ( usm_par(2,jw,iw)==0 )  THEN
+                 surf_usm_h%isroof_surf(m)    = .FALSE.
+                 surf_usm_h%thickness_wall(m) = -1.0_wp
+              ELSE
+                 surf_usm_h%isroof_surf(m)    = .TRUE.
+                 surf_usm_h%thickness_wall(m) = usm_val(2,jw,iw)
+              ENDIF
+              surf_usm_h%surface_types(m) = usm_par(5,jw,iw)
+              surf_usm_h%albedo_surf(m)   = usm_val(1,jw,iw)
+           ENDIF
+!
+!--        Find the type position
+           it = surf_usm_h%surface_types(m)
+           ip = -99999
+           DO k = 1, n_surface_types
+              IF ( surface_type_codes(k) == it )  THEN
+                 ip = k
+                 EXIT
+              ENDIF
+           ENDDO
+           IF ( ip == -99999 )  THEN
+!--           wall category not found
+              WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it,     &
+                                     ' not found  for i,j,k=', iw,jw,kw
+              CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
+           ENDIF
+!
+!--        Albedo
+           IF ( surf_usm_h%albedo_surf(m) < 0.0_wp )  THEN
+              surf_usm_h%albedo_surf(m) = surface_params(ialbedo,ip)
+           ENDIF
+!
+!--        emissivity of the wall
+           surf_usm_h%emiss_surf(m) = surface_params(iemiss,ip)
+!            
+!--        heat conductivity Î»S between air and wall ( W mâ2 Kâ1 )
+           surf_usm_h%lambda_surf(m) = surface_params(ilambdas,ip)
+!            
+!--        roughness relative to concrete
+           surf_usm_h%roughness_wall(m) = surface_params(irough,ip)
+!            
+!--        Surface skin layer heat capacity (J mâ2 Kâ1 )
+           surf_usm_h%c_surface(m) = surface_params(icsurf,ip)
+!            
+!--        wall material parameters:
+!--        thickness of the wall (m)
+!--        missing values are replaced by default value for category
+           IF ( surf_usm_h%thickness_wall(m) <= 0.001_wp )  THEN
+                surf_usm_h%thickness_wall(m) = surface_params(ithick,ip)
+           ENDIF
+!            
+!--        volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )
+           surf_usm_h%rho_c_wall(:,m) = surface_params(irhoC,ip)
+!            
+!--        thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
+           surf_usm_h%lambda_h(:,m) = surface_params(ilambdah,ip)
+
+        ENDDO
+!
+!--     For vertical surface elements ( 0 -- northward-facing, 1 -- southward-facing,
+!--     2 -- eastward-facing, 3 -- westward-facing )
+        DO  l = 0, 3
+           DO  m = 1, surf_usm_v(l)%ns
+              i  = surf_usm_v(l)%i(m)
+              j  = surf_usm_v(l)%j(m)
+              kw = surf_usm_v(l)%k(m)
+
+              IF ( l == 3 )  THEN ! westward facing
+                 iw = i
+                 jw = j
+                 ii = 6
+                 ij = 3
+              ELSEIF ( l == 2 )  THEN
+                 iw = i-1
+                 jw = j
+                 ii = 6
+                 ij = 3
+              ELSEIF ( l == 1 )  THEN
+                 iw = i
+                 jw = j
+                 ii = 12
+                 ij = 9
+              ELSEIF ( l == 0 )  THEN
+                 iw = i
+                 jw = j-1
+                 ii = 12
+                 ij = 9
+              ENDIF
+
+              IF ( kw <= usm_par(ii,jw,iw) )  THEN
+!--              pedestrant zone
+                 IF ( usm_par(ii+1,jw,iw) == 0 )  THEN
+                     surf_usm_v(l)%surface_types(m)  = pedestrant_category   !< default category for wall surface in pedestrant zone
+                     surf_usm_v(l)%albedo_surf(m)    = -1.0_wp
+                     surf_usm_v(l)%thickness_wall(m) = -1.0_wp
+                 ELSE
+                     surf_usm_v(l)%surface_types(m)  = usm_par(ii+1,jw,iw)
+                     surf_usm_v(l)%albedo_surf(m)    = usm_val(ij,jw,iw)
+                     surf_usm_v(l)%thickness_wall(m) = usm_val(ij+1,jw,iw)
+                 ENDIF
+              ELSE IF ( kw <= usm_par(ii+2,jw,iw) )  THEN
+!--              wall zone
+                 IF ( usm_par(ii+3,jw,iw) == 0 )  THEN
+                     surf_usm_v(l)%surface_types(m)  = wall_category         !< default category for wall surface
+                     surf_usm_v(l)%albedo_surf(m)    = -1.0_wp
+                     surf_usm_v(l)%thickness_wall(m) = -1.0_wp
+                 ELSE
+                     surf_usm_v(l)%surface_types(m)  = usm_par(ii+3,jw,iw)
+                     surf_usm_v(l)%albedo_surf(m)    = usm_val(ij+2,jw,iw)
+                     surf_usm_v(l)%thickness_wall(m) = usm_val(ij+3,jw,iw)
+                 ENDIF
+              ELSE IF ( kw <= usm_par(ii+4,jw,iw) )  THEN
+!--              roof zone
+                 IF ( usm_par(ii+5,jw,iw) == 0 )  THEN
+                     surf_usm_v(l)%surface_types(m)  = roof_category         !< default category for roof surface
+                     surf_usm_v(l)%albedo_surf(m)    = -1.0_wp
+                     surf_usm_v(l)%thickness_wall(m) = -1.0_wp
+                 ELSE
+                     surf_usm_v(l)%surface_types(m)  = usm_par(ii+5,jw,iw)
+                     surf_usm_v(l)%albedo_surf(m)    = usm_val(ij+4,jw,iw)
+                     surf_usm_v(l)%thickness_wall(m) = usm_val(ij+5,jw,iw)
+                 ENDIF
+              ELSE
+!--              something wrong
+                 CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 )
+              ENDIF
+
+!
+!--           Find the type position
+              it = surf_usm_v(l)%surface_types(m)
+              ip = -99999
+              DO k = 1, n_surface_types
+                 IF ( surface_type_codes(k) == it )  THEN
                     ip = k
                     EXIT
-                ENDIF
-            ENDDO
-            IF ( ip == -99999 )  THEN
-!--             wall category not found
-                WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, ' not found  for i,j,k=', iw,jw,kw
-                CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
-            ENDIF
-            
-!--         Fill out the parameters of the wall
-!--         wall surface:
-            
-!--         albedo
-            IF ( albedo_surf(l) < 0.0_wp )  THEN
-                albedo_surf(l) = surface_params(ialbedo, ip)
-            ENDIF
-            
-!--         emissivity of the wall
-            emiss_surf(l) = surface_params(iemiss, ip)
-            
-!--         heat conductivity Î»S between air and wall ( W mâ2 Kâ1 )
-            lambda_surf(l) = surface_params(ilambdas, ip)
-            
-!--         roughness relative to concrete
-            roughness_wall(l) = surface_params(irough, ip)
-            
-!--         Surface skin layer heat capacity (J mâ2 Kâ1 )
-            c_surface(l) = surface_params(icsurf, ip)
-            
-!--         wall material parameters:
-            
-!--         thickness of the wall (m)
-!--         missing values are replaced by default value for category
-            IF ( thickness_wall(l) <= 0.001_wp )  THEN
-                thickness_wall(l) = surface_params(ithick, ip)
-            ENDIF
-            
-!--         volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )
-            rho_c_wall(:,l) = surface_params(irhoC, ip)
-            
-!--         thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
-            lambda_h(:,l) = surface_params(ilambdah, ip)
-            
+                 ENDIF
+              ENDDO
+              IF ( ip == -99999 )  THEN
+!--              wall category not found
+                 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it,  &
+                                        ' not found  for i,j,k=', iw,jw,kw
+                 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
+              ENDIF
+!
+!--           Albedo
+              IF ( surf_usm_v(l)%albedo_surf(m) < 0.0_wp )  THEN
+                 surf_usm_v(l)%albedo_surf(m) = surface_params(ialbedo,ip)
+              ENDIF
+!
+!--           emissivity of the wall
+              surf_usm_v(l)%emiss_surf(m) = surface_params(iemiss,ip)
+!            
+!--           heat conductivity Î»S between air and wall ( W mâ2 Kâ1 )
+              surf_usm_v(l)%lambda_surf(m) = surface_params(ilambdas,ip)
+!            
+!--           roughness relative to concrete
+              surf_usm_v(l)%roughness_wall(m) = surface_params(irough,ip)
+!            
+!--           Surface skin layer heat capacity (J mâ2 Kâ1 )
+              surf_usm_v(l)%c_surface(m) = surface_params(icsurf,ip)
+!            
+!--           wall material parameters:
+!--           thickness of the wall (m)
+!--           missing values are replaced by default value for category
+              IF ( surf_usm_v(l)%thickness_wall(m) <= 0.001_wp )  THEN
+                   surf_usm_v(l)%thickness_wall(m) = surface_params(ithick,ip)
+              ENDIF
+!            
+!--           volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )
+              surf_usm_v(l)%rho_c_wall(:,m) = surface_params(irhoC,ip)
+!            
+!--           thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
+              surf_usm_v(l)%lambda_h(:,m) = surface_params(ilambdah,ip)
+
+           ENDDO
         ENDDO
 
@@ -3900,5 +4715,5 @@
         IMPLICIT NONE
 
-        INTEGER(iwp)                          :: i, j, k, l, d      !< running indices
+        INTEGER(iwp)                          :: i, j, k, l, d, m   !< running indices
         
         REAL(wp)                              :: pt1                !< temperature at first grid box adjacent to surface
@@ -3924,188 +4739,277 @@
         
         exn(:) = (hyp(nzub:nzut) / 100000.0_wp )**0.286_wp          !< Exner function
+!        
+!--     First, treat horizontal surface elements
+
+        DO  m = 1, surf_usm_h%ns
+!
+!--        Get indices of respective grid point
+           i = surf_usm_h%i(m)
+           j = surf_usm_h%j(m)
+           k = surf_usm_h%k(m)
+!
+!--        TODO - how to calculate lambda_surface for horizontal surfaces
+!--        (lambda_surface is set according to stratification in land surface model)
+!--        MS: ???
+           IF ( surf_usm_h%ol(m) >= 0.0_wp )  THEN
+              lambda_surface = surf_usm_h%lambda_surf(m)
+           ELSE
+              lambda_surface = surf_usm_h%lambda_surf(m)
+           ENDIF
             
-!--    
-        DO l = startenergy, endenergy
-!--         Calculate frequently used parameters
-            d = surfl(id,l)
-            k = surfl(iz,l)
-            j = surfl(iy,l)
-            i = surfl(ix,l)
-
-!--         TODO - how to calculate lambda_surface for horizontal surfaces
-!--         (lambda_surface is set according to stratification in land surface model)
-            IF ( ol(j,i) >= 0.0_wp )  THEN
-                lambda_surface = lambda_surf(l)
-            ELSE
-                lambda_surface = lambda_surf(l)
-            ENDIF
+           pt1  = pt(k,j,i)
+!
+!--        calculate rho * cp coefficient at surface layer
+           rho_cp  = cp * hyp(k) / ( r_d * pt1 * exn(k) )
+!
+!--        Calculate aerodyamic resistance. 
+!--        Calculation for horizontal surfaces follows LSM formulation
+!--        pt, us, ts are not available for the prognostic time step,
+!--        data from the last time step is used here.
+               
+           r_a = ( pt1 - t_surf_h(m) / exn(k) ) /                              &
+                 ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp )
+                
+!--        make sure that the resistance does not drop to zero
+           IF ( ABS(r_a) < 1.0E-10_wp )  r_a = 1.0E-10_wp
+                
+!--        the parameterization is developed originally for larger scales
+!--        (compare with remark in TUF-3D)
+!--        our first experiences show that the parameterization underestimates
+!--        r_a in meter resolution.
+!--        temporary solution - multiplication by magic constant :-(.
+           r_a = r_a * ra_horiz_coef
+                
+!--        factor for shf_eb
+           f_shf  = rho_cp / r_a
+        
+!--        add LW up so that it can be removed in prognostic equation
+           surf_usm_h%rad_net_l(m) = surf_usm_h%rad_in_sw(m)  -                &
+                                     surf_usm_h%rad_out_sw(m) +                &
+                                     surf_usm_h%rad_in_lw(m)  -                &
+                                     surf_usm_h%rad_out_lw(m)
+
+!--        numerator of the prognostic equation
+           coef_1 = surf_usm_h%rad_net_l(m) +                                  & 
+                   ( 3.0_wp + 1.0_wp ) * surf_usm_h%emiss_surf(m) * sigma_sb * &
+                                       t_surf_h(m) ** 4 +                      &  
+                                       f_shf * pt1 +                           &
+                                       lambda_surface * t_wall_h(nzb_wall,m)
+
+!--        denominator of the prognostic equation
+           coef_2 = 4.0_wp * surf_usm_h%emiss_surf(m) * sigma_sb *            &
+                             t_surf_h(m) ** 3                                 &
+                           + lambda_surface + f_shf / exn(k)
+
+!--        implicit solution when the surface layer has no heat capacity,
+!--        otherwise use RK3 scheme.
+           t_surf_h_p(m) = ( coef_1 * dt_3d * tsc(2) +                        &
+                             surf_usm_h%c_surface(m) * t_surf_h(m) ) /        & 
+                           ( surf_usm_h%c_surface(m) + coef_2 * dt_3d * tsc(2) ) 
+
+!--        add RK3 term
+           t_surf_h_p(m) = t_surf_h_p(m) + dt_3d * tsc(3) *                   &
+                           surf_usm_h%tt_surface_m(m)
             
-            pt1  = pt(k,j,i)
-
-!--         calculate rho * cp coefficient at surface layer
-            rho_cp  = cp * hyp(k) / ( r_d * pt1 * exn(k) )
-
-!--         calculate aerodyamic resistance. 
-            IF ( d == iroof )  THEN
-!--             calculation for horizontal surfaces follows LSM formulation
-!--             pt, us, ts are not available for the prognostic time step,
-!--             data from the last time step is used here.
+!--        calculate true tendency
+           stend = ( t_surf_h_p(m) - t_surf_h(m) - dt_3d * tsc(3) *           &
+                     surf_usm_h%tt_surface_m(m)) / ( dt_3d  * tsc(2) )
+
+!--        calculate t_surf tendencies for the next Runge-Kutta step
+           IF ( timestep_scheme(1:5) == 'runge' )  THEN
+              IF ( intermediate_timestep_count == 1 )  THEN
+                 surf_usm_h%tt_surface_m(m) = stend
+              ELSEIF ( intermediate_timestep_count <                          &
+                        intermediate_timestep_count_max )  THEN
+                 surf_usm_h%tt_surface_m(m) = -9.5625_wp * stend +            &
+                                     5.3125_wp * surf_usm_h%tt_surface_m(m)
+              ENDIF
+           ENDIF
+
+!--        in case of fast changes in the skin temperature, it is required to
+!--        update the radiative fluxes in order to keep the solution stable
+           IF ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp )  THEN
+              force_radiation_call_l = .TRUE.
+           ENDIF
+!           
+!--        for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf.
+!--        there is no equivalent surface gridpoint for vertical surfaces.
+!--        pt(k,j,i) is calculated for all directions in diffusion_s
+!--        using surface and wall heat fluxes
+           pt(k-1,j,i) = t_surf_h_p(m) / exn(k)  ! not for vertical surfaces
+
+!--        calculate fluxes
+!--        rad_net_l is never used!           
+           surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) +                &
+                                     3.0_wp * sigma_sb *                      &
+                                     t_surf_h(m)**4 - 4.0_wp * sigma_sb *     &
+                                     t_surf_h(m)**3 * t_surf_h_p(m)
+           surf_usm_h%wghf_eb(m)   = lambda_surface *                         &
+                                      ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) )
+!
+!--        ground/wall/roof surface heat flux
+           surf_usm_h%wshf_eb(m)   = - f_shf  * ( pt1 - t_surf_h_p(m) )
+!           
+!--        store kinematic surface heat fluxes for utilization in other processes
+!--        diffusion_s, surface_layer_fluxes,...
+           surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / rho_cp
+
+       ENDDO
+!
+!--    Now, treat vertical surface elements
+       DO  l = 0, 3
+          DO  m = 1, surf_usm_v(l)%ns
+!
+!--          Get indices of respective grid point
+             i = surf_usm_v(l)%i(m)
+             j = surf_usm_v(l)%j(m)
+             k = surf_usm_v(l)%k(m)
+!
+!--          TODO - how to calculate lambda_surface for horizontal (??? do you mean verical ???) surfaces
+!--          (lambda_surface is set according to stratification in land surface model).
+!--          Please note, for vertical surfaces no ol is defined, since 
+!--          stratification is not considered in this case.
+             lambda_surface = surf_usm_v(l)%lambda_surf(m)
+            
+             pt1  = pt(k,j,i)
+!
+!--          calculate rho * cp coefficient at surface layer
+             rho_cp  = cp * hyp(k) / ( r_d * pt1 * exn(k) )
+
+!--          Calculation of r_a for vertical surfaces
+!--
+!--          heat transfer coefficient for forced convection along vertical walls
+!--          follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
+!--            
+!--          H = httc (Tsfc - Tair)
+!--          httc = rw * (11.8 + 4.2 * Ueff) - 4.0
+!--            
+!--                rw: wall patch roughness relative to 1.0 for concrete
+!--                Ueff: effective wind speed
+!--                - 4.0 is a reduction of Rowley et al (1930) formulation based on
+!--                Cole and Sturrock (1977)
+!--           
+!--                Ucan: Canyon wind speed
+!--                wstar: convective velocity
+!--                Qs: surface heat flux
+!--                zH: height of the convective layer
+!--                wstar = (g/Tcan*Qs*zH)**(1./3.)
                 
-                r_a = (pt1 - t_surf(l)/exn(k)) / (ts(j,i) * us(j,i) + 1.0E-10_wp)
+!--          Effective velocity components must always 
+!--          be defined at scalar grid point. The wall normal component is 
+!--          obtained by simple linear interpolation. ( An alternative would
+!--          be an logarithmic interpolation. )
+             u1 = ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp
+             v1 = ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp
+             w1 = ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp
                 
-!--             make sure that the resistance does not drop to zero
-                IF ( ABS(r_a) < 1.0E-10_wp )  r_a = 1.0E-10_wp
-                
-!--             the parameterization is developed originally for larger scales
-!--             (compare with remark in TUF-3D)
-!--             our first experiences show that the parameterization underestimates
-!--             r_a in meter resolution.
-!--             temporary solution - multiplication by magic constant :-(.
-                r_a = r_a * ra_horiz_coef
-                
-!--             factor for shf_eb
-                f_shf  = rho_cp / r_a
-            ELSE
-!--             calculation of r_a for vertical surfaces
-!--
-!--             heat transfer coefficient for forced convection along vertical walls
-!--             follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
-!--            
-!--             H = httc (Tsfc - Tair)
-!--             httc = rw * (11.8 + 4.2 * Ueff) - 4.0
-!--            
-!--                   rw: wall patch roughness relative to 1.0 for concrete
-!--                   Ueff: effective wind speed
-!--                   - 4.0 is a reduction of Rowley et al (1930) formulation based on
-!--                   Cole and Sturrock (1977)
-!--            
-!--                   Ucan: Canyon wind speed
-!--                   wstar: convective velocity
-!--                   Qs: surface heat flux
-!--                   zH: height of the convective layer
-!--                   wstar = (g/Tcan*Qs*zH)**(1./3.)
-                
-!--             staggered grid needs to be taken into consideration
-                IF ( d == inorth )  THEN
-                    u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp
-                    v1 = v(k,j+1,i)
-                ELSE IF ( d == isouth )  THEN
-                    u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp
-                    v1 = v(k,j,i)
-                ELSE IF ( d == ieast )  THEN
-                    u1 = u(k,j,i+1)
-                    v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp
-                ELSE IF ( d == iwest )  THEN
-                    u1 = u(k,j,i)
-                    v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp
-                ELSE
-                    STOP
-                ENDIF
-                w1 = (w(k,j,i)+w(k-1,j,i))*0.5_wp
-                
-                Ueff = SQRT(u1**2 + v1**2 + w1**2)
-                httc = roughness_wall(l) * (11.8 + 4.2 * Ueff) - 4.0
-                f_shf  = httc
-            ENDIF
-        
-!--         add LW up so that it can be removed in prognostic equation
-            rad_net_l(l) = surfinsw(l) - surfoutsw(l) + surfinlw(l) - surfoutlw(l)
-
-!--         numerator of the prognostic equation
-            coef_1 = rad_net_l(l) +    &    ! coef +1 corresponds to -lwout included in calculation of radnet_l
-                     (3.0_wp+1.0_wp) * emiss_surf(l) * sigma_sb * t_surf(l) ** 4 +      &  
-                     f_shf  * pt1 +                                                     &
-                     lambda_surface * t_wall(nzb_wall,l)
-
-!--         denominator of the prognostic equation
-            coef_2 = 4.0_wp * emiss_surf(l) * sigma_sb * t_surf(l) ** 3                 &
-                         + lambda_surface + f_shf / exn(k)
-
-!--         implicit solution when the surface layer has no heat capacity,
-!--         otherwise use RK3 scheme.
-            t_surf_p(l) = ( coef_1 * dt_3d * tsc(2) + c_surface(l) * t_surf(l) ) /      & 
-                              ( c_surface(l) + coef_2 * dt_3d * tsc(2) ) 
-
-!--         add RK3 term
-            t_surf_p(l) = t_surf_p(l) + dt_3d * tsc(3) * tt_surface_m(l)
+             Ueff   = SQRT( u1**2 + v1**2 + w1**2 )
+             httc   = surf_usm_v(l)%roughness_wall(m) *                       &
+                       ( 11.8_wp + 4.2_wp * Ueff ) - 4.0_wp
+             f_shf  = httc
+
+!--          add LW up so that it can be removed in prognostic equation
+             surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_in_sw(m)  -        &
+                                          surf_usm_v(l)%rad_out_sw(m) +        &
+                                          surf_usm_v(l)%rad_in_lw(m)  -        &
+                                          surf_usm_v(l)%rad_out_lw(m)
+
+!--           numerator of the prognostic equation
+              coef_1 = surf_usm_v(l)%rad_net_l(m) +                            & ! coef +1 corresponds to -lwout included in calculation of radnet_l
+               ( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emiss_surf(m) * sigma_sb *  &
+                                     t_surf_v(l)%t(m) ** 4 +                   &  
+                                     f_shf * pt1 +                             &
+                                     lambda_surface * t_wall_v(l)%t(nzb_wall,m)
+
+!--           denominator of the prognostic equation
+              coef_2 = 4.0_wp * surf_usm_v(l)%emiss_surf(m) * sigma_sb *       &
+                                t_surf_v(l)%t(m) ** 3                          &
+                              + lambda_surface + f_shf / exn(k)
+
+!--           implicit solution when the surface layer has no heat capacity,
+!--           otherwise use RK3 scheme.
+              t_surf_v_p(l)%t(m) = ( coef_1 * dt_3d * tsc(2) +                 &
+                             surf_usm_v(l)%c_surface(m) * t_surf_v(l)%t(m) ) / & 
+                           ( surf_usm_v(l)%c_surface(m) + coef_2 * dt_3d * tsc(2) ) 
+
+!--           add RK3 term
+              t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) *       &
+                                surf_usm_v(l)%tt_surface_m(m)
             
-!--         calculate true tendency
-            stend = (t_surf_p(l) - t_surf(l) - dt_3d * tsc(3) * tt_surface_m(l)) / (dt_3d  * tsc(2))
-
-!--         calculate t_surf tendencies for the next Runge-Kutta step
-            IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
-                    tt_surface_m(l) = stend
-                ELSEIF ( intermediate_timestep_count <                                  &
-                         intermediate_timestep_count_max )  THEN
-                    tt_surface_m(l) = -9.5625_wp * stend + 5.3125_wp                    &
-                                       * tt_surface_m(l)
-                ENDIF
-            ENDIF
-
-!--         in case of fast changes in the skin temperature, it is required to
-!--         update the radiative fluxes in order to keep the solution stable
-            IF ( ABS( t_surf_p(l) - t_surf(l) ) > 1.0_wp )  THEN
-               force_radiation_call_l = .TRUE.
-            ENDIF
-            
-!--         for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf.
-!--         there is no equivalent surface gridpoint for vertical surfaces.
-!--         pt(k,j,i) is calculated for all directions in diffusion_s
-!--         using surface and wall heat fluxes
-            IF ( d == iroof )  THEN
-               pt(nzb_s_inner(j,i),j,i) = t_surf_p(l) / exn(k)
-            ENDIF
-
-!--         calculate fluxes
-!--         rad_net_l is never used!           
-            rad_net_l(l)     = rad_net_l(l) + 3.0_wp * sigma_sb                         &
-                                * t_surf(l)**4 - 4.0_wp * sigma_sb                      &
-                                * t_surf(l)**3 * t_surf_p(l)
-            wghf_eb(l)       = lambda_surface * (t_surf_p(l) - t_wall(nzb_wall,l))
-
-!--         ground/wall/roof surface heat flux
-            wshf_eb(l)  = - f_shf  * ( pt1 - t_surf_p(l) )
-            
-!--         store kinematic surface heat fluxes for utilization in other processes
-!--         diffusion_s, surface_layer_fluxes,...
-            IF ( d == iroof )  THEN
-!--             shf is used in diffusion_s and also
-!--             for calculation of surface layer fluxes
-!--             update for horizontal surfaces
-                shf(j,i) = wshf_eb(l) / rho_cp
-            ELSE
-!--             surface heat flux for vertical surfaces
-!--             used in diffusion_s
-                wshf(l) = wshf_eb(l) / rho_cp
-            ENDIF
+!--           calculate true tendency
+              stend = ( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) - dt_3d * tsc(3) *&
+                        surf_usm_v(l)%tt_surface_m(m) ) / ( dt_3d  * tsc(2) )
+
+!--           calculate t_surf tendencies for the next Runge-Kutta step
+              IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                 IF ( intermediate_timestep_count == 1 )  THEN
+                    surf_usm_v(l)%tt_surface_m(m) = stend
+                 ELSEIF ( intermediate_timestep_count <                        &
+                          intermediate_timestep_count_max )  THEN
+                    surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend +       &
+                                     5.3125_wp * surf_usm_h%tt_surface_m(m)
+                 ENDIF
+              ENDIF
+
+!--           in case of fast changes in the skin temperature, it is required to
+!--           update the radiative fluxes in order to keep the solution stable
+              IF ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp )  THEN
+                 force_radiation_call_l = .TRUE.
+              ENDIF
+
+!--           calculate fluxes
+!--           rad_net_l is never used!           
+              surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_net_l(m) +        &
+                                     3.0_wp * sigma_sb *                       &
+                                     t_surf_v(l)%t(m)**4 - 4.0_wp * sigma_sb * &
+                                     t_surf_v(l)%t(m)**3 * t_surf_v_p(l)%t(m)
+
+              surf_usm_v(l)%wghf_eb(m)   = lambda_surface *                    &
+                                     ( t_surf_v_p(l)%t(m) - t_wall_v(l)%t(nzb_wall,m) )
+
+!--           ground/wall/roof surface heat flux
+              surf_usm_v(l)%wshf_eb(m)   = - f_shf  * ( pt1 - t_surf_v_p(l)%t(m) )
+
+!           
+!--           store kinematic surface heat fluxes for utilization in other processes
+!--           diffusion_s, surface_layer_fluxes,...
+              surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / rho_cp
+
+           ENDDO
 
         ENDDO
-        
-        
+!
+!--     Add-up anthropogenic heat, for now only at upward-facing surfaces
         IF ( usm_anthropogenic_heat  .AND.  &
              intermediate_timestep_count == intermediate_timestep_count_max )  THEN
-!--         application of the additional anthropogenic heat sources
-!--         we considere the traffic for now so all heat is absorbed
-!--         to the first layer, generalization would be worth
+!--        application of the additional anthropogenic heat sources
+!--        we considere the traffic for now so all heat is absorbed
+!--        to the first layer, generalization would be worth. 
             
-!--         calculation of actual profile coefficient
-!--         ??? check time_since_reference_point ???
-            dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
-            dhour = INT(dtime/3600.0_wp)
-!--         linear interpolation of coeficient
-            acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1)
-            DO i = nxl, nxr
-                DO j = nys, nyn
-                    IF ( aheat(j,i) > 0.0_wp )  THEN
-!--                     TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d 
-!--                     given to anthropogenic heat aheat*acoef (W*m-2)
-!--                     k = nzb_s_inner(j,i)+1
-!--                     pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz)
-!--                     Instead of this, we can adjust shf in case AH only at surface
-                        shf(j,i) = shf(j,i) + aheat(j,i)*acoef * ddx * ddy / rho_cp
-                    ENDIF
-                ENDDO
-            ENDDO
+!--        calculation of actual profile coefficient
+!--        ??? check time_since_reference_point ???
+           dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
+           dhour = INT(dtime/3600.0_wp)
+!--        linear interpolation of coeficient
+           acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1)
+
+           DO m = 1, surf_usm_h%ns
+!
+!--           Get indices of respective grid point
+              i = surf_usm_h%i(m)
+              j = surf_usm_h%j(m)
+              k = surf_usm_h%k(m)
+
+              IF ( aheat(j,i) > 0.0_wp )  THEN
+!--              TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d 
+!--              given to anthropogenic heat aheat*acoef (W*m-2)
+!--              k = nzb_s_inner(j,i)+1
+!--              pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz)
+!--              Instead of this, we can adjust shf in case AH only at surface
+                 surf_usm_h%shf(m) = surf_usm_h%shf(m) +                       &
+                                   aheat(j,i) * acoef * ddx * ddy / rho_cp
+              ENDIF
+           ENDDO
+
         ENDIF
         
@@ -4113,16 +5017,15 @@
 !--     get the borders from neighbours
         CALL exchange_horiz( pt, nbgp )
-        CALL exchange_horiz_2d( shf )
-
-
-!--    calculation of force_radiation_call:
-!--    Make logical OR for all processes.
-!--    Force radiation call if at least one processor forces it.
-       IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )          &
-       THEN
+
+
+!--     calculation of force_radiation_call:
+!--     Make logical OR for all processes.
+!--     Force radiation call if at least one processor forces it.
+        IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )&
+        THEN
 #if defined( __parallel )
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-              CALL mpi_allreduce( force_radiation_call_l, force_radiation_call,         &
-                                  1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
+          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
+                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
 #else
           force_radiation_call = force_radiation_call_l
@@ -4130,5 +5033,4 @@
           force_radiation_call_l = .FALSE.
        ENDIF
-
 
     END SUBROUTINE usm_surface_energy_balance
@@ -4149,77 +5051,32 @@
       
 #if defined( __nopointer )
-       t_surf    = t_surf_p
-       t_wall    = t_wall_p
+       t_surf_h    = t_surf_h_p
+       t_wall_h    = t_wall_h_p
+       t_surf_v    = t_surf_v_p
+       t_wall_v    = t_wall_v_p
 #else
        SELECT CASE ( mod_count )
           CASE ( 0 )
-             t_surf  => t_surf_1; t_surf_p  => t_surf_2
-             t_wall     => t_wall_1;    t_wall_p     => t_wall_2
+!
+!--          Horizontal surfaces
+             t_surf_h  => t_surf_h_1; t_surf_h_p  => t_surf_h_2
+             t_wall_h     => t_wall_h_1;    t_wall_h_p     => t_wall_h_2
+!
+!--          Vertical surfaces
+             t_surf_v  => t_surf_v_1; t_surf_v_p  => t_surf_v_2
+             t_wall_v     => t_wall_v_1;    t_wall_v_p     => t_wall_v_2
           CASE ( 1 )
-             t_surf  => t_surf_2; t_surf_p  => t_surf_1
-             t_wall     => t_wall_2;    t_wall_p     => t_wall_1
+!
+!--          Horizontal surfaces
+             t_surf_h  => t_surf_h_2; t_surf_h_p  => t_surf_h_1
+             t_wall_h     => t_wall_h_2;    t_wall_h_p     => t_wall_h_1
+!
+!--          Vertical surfaces
+             t_surf_v  => t_surf_v_2; t_surf_v_p  => t_surf_v_1
+             t_wall_v     => t_wall_v_2;    t_wall_v_p     => t_wall_v_1
        END SELECT
 #endif
         
     END SUBROUTINE usm_swap_timelevel
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!
-!> This function applies the kinematic wall heat fluxes
-!> for walls in four directions for all gridboxes in urban layer.
-!> It is called out from subroutine prognostic_equations.
-!> TODO Compare performance with cycle runnig l=startwall,endwall...
-!------------------------------------------------------------------------------!
-    SUBROUTINE usm_wall_heat_flux
-    
-        IMPLICIT NONE
-
-        INTEGER(iwp)              ::  i,j,k,d,l             !< running indices
-        
-        DO l = startenergy, endenergy
-            j = surfl(iy,l)
-            i = surfl(ix,l)
-            k = surfl(iz,l)
-            d = surfl(id,l)
-            tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)
-        ENDDO
-
-    END SUBROUTINE usm_wall_heat_flux
- 
- 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!
-!> This function applies the kinematic wall heat fluxes
-!> for walls in four directions around the gridbox i,j.
-!> It is called out from subroutine prognostic_equations.
-!------------------------------------------------------------------------------!
-    SUBROUTINE usm_wall_heat_flux_ij(i,j) 
-    
-        IMPLICIT NONE
-
-        INTEGER(iwp), INTENT(in)  ::  i,j                   !< indices of grid box
-        INTEGER(iwp)              ::  ii,jj,k,d,l
-        
-        DO l = startenergy, endenergy
-            jj = surfl(iy,l)
-            ii = surfl(ix,l)
-            IF ( ii == i  .AND.  jj == j ) THEN
-               k = surfl(iz,l)
-               IF ( k >=  nzb_s_inner(j,i)+1  .AND.  k <=  nzb_s_outer(j,i) ) THEN
-                  d = surfl(id,l)
-                  IF ( d >= 1 .and. d <= 4 )   THEN
-                     tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)
-                  ENDIF
-               ENDIF
-            ENDIF
-        ENDDO
-
-    END SUBROUTINE usm_wall_heat_flux_ij
- 
 
 !------------------------------------------------------------------------------!
@@ -4242,15 +5099,64 @@
        DO  i = 0, io_blocks-1
           IF ( i == io_group )  THEN
-             WRITE ( 14 )  't_surf                        '
+
+             WRITE ( 14 )  't_surf_h                      '
 #if defined( __nopointer )             
-             WRITE ( 14 )  t_surf
+             WRITE ( 14 )  t_surf_h
 #else
-             WRITE ( 14 )  t_surf_1
+             WRITE ( 14 )  t_surf_h_1
 #endif
-             WRITE ( 14 )  't_wall                        '
+             WRITE ( 14 )  't_surf_v(0)                   '
 #if defined( __nopointer )             
-             WRITE ( 14 )  t_wall
+             WRITE ( 14 )  t_surf_v(0)%t
 #else
-             WRITE ( 14 )  t_wall_1
+             WRITE ( 14 )  t_surf_v_1(0)%t
+#endif
+             WRITE ( 14 )  't_surf_v(1)                   '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_surf_v(1)%t
+#else
+             WRITE ( 14 )  t_surf_v_1(1)%t
+#endif
+             WRITE ( 14 )  't_surf_v(2)                   '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_surf_v(2)%t
+#else
+             WRITE ( 14 )  t_surf_v_1(2)%t
+#endif
+             WRITE ( 14 )  't_surf_v(3)                   '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_surf_v(3)%t
+#else
+             WRITE ( 14 )  t_surf_v_1(3)%t
+#endif
+             WRITE ( 14 )  't_wall_h                      '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_wall_h
+#else
+             WRITE ( 14 )  t_wall_h_1
+#endif
+             WRITE ( 14 )  't_wall_v(0)                   '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_wall_v(0)%t
+#else
+             WRITE ( 14 )  t_wall_v_1(0)%t
+#endif
+             WRITE ( 14 )  't_wall_v(1)                   '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_wall_v(1)%t
+#else
+             WRITE ( 14 )  t_wall_v_1(1)%t
+#endif
+             WRITE ( 14 )  't_wall_v(2)                   '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_wall_v(2)%t
+#else
+             WRITE ( 14 )  t_wall_v_1(2)%t
+#endif
+             WRITE ( 14 )  't_wall_v(3)                   '
+#if defined( __nopointer )             
+             WRITE ( 14 )  t_wall_v(3)%t
+#else
+             WRITE ( 14 )  t_wall_v_1(3)%t
 #endif
              WRITE ( 14 )  '*** end usm ***               '
