Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3274)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3337)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw*
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -349,5 +352,5 @@
                spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
                average_count_3d, varnamelength, urban_surface,                 &
-               plant_canopy
+               plant_canopy, dz
 
     USE cpulog,                                                                &
@@ -381,11 +384,11 @@
                surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
                rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
-               surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw,             &
-               iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,             &
+               surfl, nsurfl, pcbinsw, pcbinlw, pcbinswdir,                    &
+               pcbinswdif, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, &
                inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
-               iz, iy, ix, idir, jdir, kdir,  nsurf_type, nsurf, idsvf, ndsvf, &
-               iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,           &
+               iz, iy, ix,  nsurf, idsvf, ndsvf,                               &
                idcsf, ndcsf, kdcsf, pct,                                       &
-               startland, endland, startwall, endwall, skyvf, skyvft
+               startland, endland, startwall, endwall, skyvf, skyvft, nzub,    &
+               nzut, nzpt, npcbl, pcbl
 
     USE statistics,                                                            &
@@ -414,4 +417,5 @@
     INTEGER(iwp) ::  pedestrian_category = 2         !< default category for wall surface in pedestrian zone
     INTEGER(iwp) ::  roof_category = 2               !< default category for root surface
+    REAL(wp)     ::  roughness_concrete = 0.001_wp   !< roughness length of average concrete surface
 !
 !-- Indices of input attributes for (above) ground floor level 
@@ -605,4 +609,9 @@
     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            ::  pcbinlw_av       !< Average of pcbinlw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
     
 
@@ -1260,9 +1269,10 @@
 !--     find the real name of the variable
         ids = -1
+        l = -1
         var = TRIM(variable)
         DO i = 0, nd-1
             k = len(TRIM(var))
             j = len(TRIM(dirname(i)))
-            IF ( var(k-j+1:k) == dirname(i) )  THEN
+            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
                 ids = i
                 idsint = dirint(ids)
@@ -1271,4 +1281,8 @@
             ENDIF
         ENDDO
+        l = idsint - 2  ! horisontal direction index - terible hack !
+        IF ( l < 0 .OR. l > 3 ) THEN
+           l = -1
+        END IF
         IF ( ids == -1 )  THEN
             var = TRIM(variable)
@@ -1311,14 +1325,13 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
+                    IF ( l == -1 .AND. .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
+                    ELSE
+                        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
                     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' )
@@ -1398,30 +1411,63 @@
                     ENDIF
                                     
+                CASE ( 'usm_rad_pc_inlw' )
+!--                 array of of lw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinlw_av) )  THEN
+                        ALLOCATE( pcbinlw_av(1:npcbl) )
+                        pcbinlw_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_insw' )
+!--                 array of of sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinsw_av) )  THEN
+                        ALLOCATE( pcbinsw_av(1:npcbl) )
+                        pcbinsw_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_inswdir' )
+!--                 array of of direct sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinswdir_av) )  THEN
+                        ALLOCATE( pcbinswdir_av(1:npcbl) )
+                        pcbinswdir_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_inswdif' )
+!--                 array of of diffuse sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinswdif_av) )  THEN
+                        ALLOCATE( pcbinswdif_av(1:npcbl) )
+                        pcbinswdif_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_inswref' )
+!--                 array of of reflected sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinswref_av) )  THEN
+                        ALLOCATE( pcbinswref_av(1:npcbl) )
+                        pcbinswref_av = 0.0_wp
+                    ENDIF
+                                   
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
+                    IF ( l == -1 .AND. .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
+                    ELSE
                        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
+                    ENDIF
 
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces
 !--                 land surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
+                    IF ( l == -1 .AND. .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
+                    ELSE
                        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
+                    ENDIF
 !
 !--             Please note, the following output quantities belongs to the 
@@ -1431,157 +1477,145 @@
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
+                    IF ( l == -1 .AND. .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
+                    ELSE
                        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
+                    ENDIF
 
                 CASE ( 'usm_wghf_window' )
 !--                 array of heat flux from window ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
                         ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
                         surf_usm_h%wghf_eb_window_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wghf_green' )
 !--                 array of heat flux from green ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
                         ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
                         surf_usm_h%wghf_eb_green_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_iwghf' )
 !--                 array of heat flux from indoor ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
                         ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
                         surf_usm_h%iwghf_eb_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%iwghf_eb_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_iwghf_window' )
 !--                 array of heat flux from indoor window ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
                         ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
                         surf_usm_h%iwghf_eb_window_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
+                    IF ( l == -1 .AND. .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
+                    ELSE
                        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
+                    ENDIF
 
                 CASE ( 'usm_t_surf_window' )
 !--                 surface temperature for window surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
                         ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
                         surf_usm_h%t_surf_window_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_surf_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_green' )
 !--                 surface temperature for green surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
                         ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
                         surf_usm_h%t_surf_green_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_green_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_surf_green_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
                 
                 CASE ( 'usm_t_surf_10cm' )
 !--                 near surface temperature for whole surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
                         ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) )
                         surf_usm_h%t_surf_10cm_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_surf_10cm_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for iwl layer of walls and land
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
+                    IF ( l == -1 .AND. .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
+                    ELSE
                        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
+                    ENDIF
 
                 CASE ( 'usm_t_window' )
 !--                 window temperature for iwl layer of walls and land
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
                         ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
                         surf_usm_h%t_window_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_green' )
 !--                 green temperature for iwl layer of walls and land
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
                         ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
                         surf_usm_h%t_green_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_green_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_green_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                CASE DEFAULT
@@ -1596,10 +1630,11 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%rad_net_av(m) =                        &
@@ -1607,6 +1642,6 @@
                                           surf_usm_v(l)%rad_net_l(m)
                        ENDDO
-                    ENDDO
-                    
+                    ENDIF
+
                 CASE ( 'usm_rad_insw' )
 !--                 array of sw radiation falling to surface after i-th reflection
@@ -1700,12 +1735,30 @@
                     ENDDO
                     
+                CASE ( 'usm_rad_pc_inlw' )
+                    pcbinlw_av(:) = pcbinlw_av(:) + pcbinlw(:)
+                    
+                CASE ( 'usm_rad_pc_insw' )
+                    pcbinsw_av(:) = pcbinsw_av(:) + pcbinsw(:)
+                    
+                CASE ( 'usm_rad_pc_inswdir' )
+                    pcbinswdir_av(:) = pcbinswdir_av(:) + pcbinswdir(:)
+                    
+                CASE ( 'usm_rad_pc_inswdif' )
+                    pcbinswdif_av(:) = pcbinswdif_av(:) + pcbinswdif(:)
+                    
+                CASE ( 'usm_rad_pc_inswref' )
+                    pcbinswref_av(:) = pcbinswref_av(:) + pcbinsw(:)     &
+                                                        - pcbinswdir(:)  &
+                                                        - pcbinswdif(:)
+                    
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%surfhf_av(m) =                         &
@@ -1713,14 +1766,15 @@
                                           surf_usm_v(l)%surfhf(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wshf_eb_av(m) =                        &
@@ -1728,14 +1782,15 @@
                                           surf_usm_v(l)%wshf_eb(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_av(m) =                        &
@@ -1743,14 +1798,15 @@
                                           surf_usm_v(l)%wghf_eb(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf_window' )
 !--                 array of heat flux from window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_window_av(m) =                              &
-                                          surf_usm_h%wghf_eb_window_av(m) +           &
-                                          surf_usm_h%wghf_eb_window(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_window_av(m) =                              &
+                                             surf_usm_h%wghf_eb_window_av(m) +           &
+                                             surf_usm_h%wghf_eb_window(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_window_av(m) =                        &
@@ -1758,14 +1814,15 @@
                                           surf_usm_v(l)%wghf_eb_window(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wghf_green' )
 !--                 array of heat flux from green ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_green_av(m) =                              &
-                                          surf_usm_h%wghf_eb_green_av(m) +           &
-                                          surf_usm_h%wghf_eb_green(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_green_av(m) =                              &
+                                             surf_usm_h%wghf_eb_green_av(m) +           &
+                                             surf_usm_h%wghf_eb_green(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_green_av(m) =                        &
@@ -1773,14 +1830,15 @@
                                           surf_usm_v(l)%wghf_eb_green(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_iwghf' )
 !--                 array of heat flux from indoor ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_av(m) +           &
-                                          surf_usm_h%iwghf_eb(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_av(m) +           &
+                                             surf_usm_h%iwghf_eb(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_av(m) =                        &
@@ -1788,14 +1846,15 @@
                                           surf_usm_v(l)%iwghf_eb(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_iwghf_window' )
 !--                 array of heat flux from indoor window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_window_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_window_av(m) +           &
-                                          surf_usm_h%iwghf_eb_window(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_window_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_window_av(m) +           &
+                                             surf_usm_h%iwghf_eb_window(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
@@ -1803,14 +1862,15 @@
                                           surf_usm_v(l)%iwghf_eb_window(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_av(m) =                         &
@@ -1818,14 +1878,15 @@
                                           t_surf_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_window' )
 !--                 surface temperature for window surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_window_av(m) =                               & 
-                                          surf_usm_h%t_surf_window_av(m) +            &
-                                          t_surf_window_h(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_window_av(m) =                               &
+                                             surf_usm_h%t_surf_window_av(m) +            &
+                                             t_surf_window_h(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_window_av(m) =                         &
@@ -1833,14 +1894,15 @@
                                           t_surf_window_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_green' )
 !--                 surface temperature for green surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_green_av(m) =                               & 
-                                          surf_usm_h%t_surf_green_av(m) +            &
-                                          t_surf_green_h(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_green_av(m) =                               &
+                                             surf_usm_h%t_surf_green_av(m) +            &
+                                             t_surf_green_h(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_green_av(m) =                         &
@@ -1848,14 +1910,15 @@
                                           t_surf_green_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                 
                 CASE ( 'usm_t_surf_10cm' )
 !--                 near surface temperature for whole surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_10cm_av(m) =                               & 
-                                          surf_usm_h%t_surf_10cm_av(m) +            &
-                                          t_surf_10cm_h(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_10cm_av(m) =                               &
+                                             surf_usm_h%t_surf_10cm_av(m) +            &
+                                             t_surf_10cm_h(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_10cm_av(m) =                         &
@@ -1863,15 +1926,16 @@
                                           t_surf_10cm_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                     
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for  iwl layer of walls and land
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_wall_av(iwl,m) =                     &
@@ -1879,14 +1943,15 @@
                                           t_wall_v(l)%t(iwl,m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_window' )
 !--                 window temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_window_av(iwl,m) =                           &
-                                          surf_usm_h%t_window_av(iwl,m) +        &
-                                          t_window_h(iwl,m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_window_av(iwl,m) =                           &
+                                             surf_usm_h%t_window_av(iwl,m) +        &
+                                             t_window_h(iwl,m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_window_av(iwl,m) =                     &
@@ -1894,14 +1959,15 @@
                                           t_window_v(l)%t(iwl,m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_green' )
 !--                 green temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_green_av(iwl,m) =                           &
-                                          surf_usm_h%t_green_av(iwl,m) +        &
-                                          t_green_h(iwl,m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_green_av(iwl,m) =                           &
+                                             surf_usm_h%t_green_av(iwl,m) +        &
+                                             t_green_h(iwl,m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_green_av(iwl,m) =                     &
@@ -1909,5 +1975,5 @@
                                           t_green_v(l)%t(iwl,m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE DEFAULT
@@ -1922,10 +1988,11 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%rad_net_av(m) =                        &
@@ -1933,5 +2000,5 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_rad_insw' )
@@ -2023,12 +2090,28 @@
                     ENDDO
                     
+                CASE ( 'usm_rad_pc_inlw' )
+                    pcbinlw_av(:) = pcbinlw_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_insw' )
+                    pcbinsw_av(:) = pcbinsw_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_inswdir' )
+                    pcbinswdir_av(:) = pcbinswdir_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_inswdif' )
+                    pcbinswdif_av(:) = pcbinswdif_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_inswref' )
+                    pcbinswref_av(:) = pcbinswref_av(:) / REAL( average_count_3d, kind=wp )
+                    
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%surfhf_av(m) =                         &
@@ -2036,14 +2119,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wshf_eb_av(m) =                        &
@@ -2051,14 +2135,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_av(m) =                        &
@@ -2066,14 +2151,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf_window' )
 !--                 array of heat flux from window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_window_av(m) =                              &
-                                          surf_usm_h%wghf_eb_window_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_window_av(m) =                              &
+                                             surf_usm_h%wghf_eb_window_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_window_av(m) =                        &
@@ -2081,14 +2167,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wghf_green' )
 !--                 array of heat flux from green ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_green_av(m) =                              &
-                                          surf_usm_h%wghf_eb_green_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_green_av(m) =                              &
+                                             surf_usm_h%wghf_eb_green_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_green_av(m) =                        &
@@ -2096,14 +2183,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_iwghf' )
 !--                 array of heat flux from indoor ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_av(m) =                        &
@@ -2111,14 +2199,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_iwghf_window' )
 !--                 array of heat flux from indoor window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_window_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_window_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_window_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_window_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
@@ -2126,14 +2215,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_av(m) =                         &
@@ -2141,14 +2231,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_window' )
 !--                 surface temperature for window surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_window_av(m) =                               & 
-                                          surf_usm_h%t_surf_window_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_window_av(m) =                               &
+                                             surf_usm_h%t_surf_window_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_window_av(m) =                         &
@@ -2156,14 +2247,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_green' )
 !--                 surface temperature for green surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_green_av(m) =                               & 
-                                          surf_usm_h%t_surf_green_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_green_av(m) =                               &
+                                             surf_usm_h%t_surf_green_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_green_av(m) =                         &
@@ -2171,14 +2263,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_10cm' )
 !--                 near surface temperature for whole surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_10cm_av(m) =                               & 
-                                          surf_usm_h%t_surf_10cm_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_10cm_av(m) =                               &
+                                             surf_usm_h%t_surf_10cm_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_10cm_av(m) =                         &
@@ -2186,14 +2279,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for  iwl layer of walls and land
-                    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
+                    IF ( l == -1 ) THEN
+                       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
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_wall_av(iwl,m) =                     &
@@ -2201,14 +2295,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_window' )
 !--                 window temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_window_av(iwl,m) =                           &
-                                          surf_usm_h%t_window_av(iwl,m) /        &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_window_av(iwl,m) =                           &
+                                             surf_usm_h%t_window_av(iwl,m) /        &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_window_av(iwl,m) =                     &
@@ -2216,14 +2311,15 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_green' )
 !--                 green temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_green_av(iwl,m) =                           &
-                                          surf_usm_h%t_green_av(iwl,m) /        &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_green_av(iwl,m) =                           &
+                                             surf_usm_h%t_green_av(iwl,m) /        &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_green_av(iwl,m) =                     &
@@ -2231,5 +2327,5 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
 
@@ -2309,12 +2405,17 @@
              var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
              var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR.     &
-             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' )  THEN
+             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR.         &
+             var(1:17) == 'usm_surfwintrans_' )  THEN
             unit = 'W/m2'
-        ELSE IF ( var(1:10) == 'usm_t_surf'   .OR.  var(1:10) == 'usm_t_wall' .OR.         &
+        ELSE IF ( var(1:10) == 'usm_t_surf'   .OR.  var(1:10) == 'usm_t_wall' .OR.        &
                   var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR.  &
                   var(1:16) == 'usm_t_surf_green'  .OR.                                   &
                   var(1:11) == 'usm_t_green' .OR.                                         &
-                  var(1:15) == 'usm_t_surf_10cm')  THEN
+                  var(1:15) == 'usm_t_surf_10cm' )  THEN
             unit = 'K'
+        ELSE IF ( var == 'usm_rad_pc_inlw'  .OR.  var == 'usm_rad_pc_insw'  .OR.          &
+                  var == 'usm_rad_pc_inswdir'  .OR.  var == 'usm_rad_pc_inswdif'  .OR.    &
+                  var == 'usm_rad_pc_inswref' )  THEN
+            unit = 'W'
         ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              & 
                   var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
@@ -2417,5 +2518,5 @@
         INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirstart
         INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
-        INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt
+        INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt,ipcgb
         INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
 
@@ -2431,5 +2532,5 @@
             k = len(TRIM(var))
             j = len(TRIM(dirname(i)))
-            IF ( var(k-j+1:k) == dirname(i) )  THEN
+            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
                 ids = i
                 idsint = dirint(ids)
@@ -2802,4 +2903,56 @@
                    ENDIF
                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inlw' )
+!--           array of lw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_insw' )
+!--           array of sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inswdir' )
+!--           array of direct sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inswdif' )
+!--           array of diffuse sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inswref' )
+!--           array of reflected sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)      &
+                                                                              - pcbinswdir(ipcgb) &
+                                                                              - pcbinswdif(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswref_av(ipcgb)
+                  ENDIF
               ENDDO
  
@@ -3356,5 +3509,5 @@
           CASE DEFAULT
               found = .FALSE.
-              
+              RETURN
         END SELECT
 
@@ -3399,5 +3552,7 @@
              var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.       &
              var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.       &
-             var(1:11) == 'usm_rad_hf_'  .OR.                                               &
+             var(1:11) == 'usm_rad_hf_'  .OR.  var == 'usm_rad_pc_inlw'  .OR.               &
+             var == 'usm_rad_pc_insw'  .OR.  var == 'usm_rad_pc_inswdir'  .OR.              &
+             var == 'usm_rad_pc_inswdif'  .OR.  var == 'usm_rad_pc_inswref'  .OR.           &
              var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
              var(1:16) == 'usm_wghf_window_'  .OR. var(1:15) == 'usm_wghf_green_' .OR.      &
@@ -4965,5 +5120,5 @@
                                  * wintend(nzb_wall:nzt_wall) + tsc(3)            &
                                  * surf_usm_h%tt_window_m(nzb_wall:nzt_wall,m) )   
-           
+
 !
 !--        calculate t_wall tendencies for the next Runge-Kutta step
@@ -5168,4 +5323,10 @@
               ENDIF
            ENDDO
+!!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!
+!           t_window_v_p(l)%t = t_wall_v_p(l)%t
+!           surf_usm_v(l)%tt_window_m  = surf_usm_v(l)%tt_wall_m
+!           t_green_v_p(l)%t = t_wall_v_p(l)%t
+!           surf_usm_v(l)%tt_green_m  = surf_usm_v(l)%tt_wall_m
+!!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!
         ENDDO
 
@@ -5320,4 +5481,5 @@
                            naheatlayers,                                       &
                            pedestrian_category,                                &
+                           roughness_concrete,                                 &
                            read_wall_temp_3d,                                  &
                            roof_category,                                      &
@@ -5337,4 +5499,5 @@
                            naheatlayers,                                       &
                            pedestrian_category,                                &
+                           roughness_concrete,                                 &
                            read_wall_temp_3d,                                  &
                            roof_category,                                      &
@@ -6937,8 +7100,25 @@
            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 )
+!--           land/roof category not found
+              WRITE (9,"(A,I5,A,3I5)") 'land/roof category ', it,     &
+                                       ' not found  for i,j,k=', iw,jw,kw
+              FLUSH(9)
+              IF ( surf_usm_h%isroof_surf(m) ) THEN
+                 category = roof_category
+              ELSE
+                 category = land_category
+              ENDIF
+              DO k = 1, n_surface_types
+                 IF ( surface_type_codes(k) == roof_category ) THEN
+                    ip = k
+                    EXIT
+                 ENDIF
+              ENDDO
+              IF ( ip == -99999 )  THEN
+!--              default land/roof category not found
+                 WRITE (9,"(A,I5,A,3I5)") 'Default land/roof category', category, ' not found!'
+                 FLUSH(9)
+                 ip = 1
+              ENDIF
            ENDIF
 !
@@ -7036,4 +7216,7 @@
                  surf_usm_v(l)%albedo(:,m)    = -1.0_wp
                  surf_usm_v(l)%thickness_wall(m) = -1.0_wp
+                 surf_usm_v(l)%thickness_window(m)   = -1.0_wp
+                 surf_usm_v(l)%thickness_green(m)    = -1.0_wp
+                 surf_usm_v(l)%transmissivity(m)  = -1.0_wp
               ELSE IF ( kw <= usm_par(ii,jw,iw) )  THEN
 !--                 pedestrian zone
@@ -7089,4 +7272,12 @@
               ELSE
 !
+                 WRITE(9,*) 'Problem reading USM data:'
+                 WRITE(9,*) l,i,j,kw,get_topography_top_index_ji( j, i, 's' )
+                 WRITE(9,*) ii,iw,jw,kw,get_topography_top_index_ji( jw, iw, 's' )
+                 WRITE(9,*) usm_par(ii,jw,iw),usm_par(ii+1,jw,iw)
+                 WRITE(9,*) usm_par(ii+2,jw,iw),usm_par(ii+3,jw,iw)
+                 WRITE(9,*) usm_par(ii+4,jw,iw),usm_par(ii+5,jw,iw)
+                 WRITE(9,*) kw,roof_height_limit,wall_category,roof_category
+                 FLUSH(9)
 !--              supply the default category
                  IF ( kw <= roof_height_limit ) THEN
@@ -7097,4 +7288,7 @@
                  surf_usm_v(l)%albedo(:,m)    = -1.0_wp
                  surf_usm_v(l)%thickness_wall(m) = -1.0_wp
+                 surf_usm_v(l)%thickness_window(m) = -1.0_wp
+                 surf_usm_v(l)%thickness_green(m) = -1.0_wp
+                 surf_usm_v(l)%transmissivity(m)  = -1.0_wp
               ENDIF
 !
@@ -7110,8 +7304,22 @@
               IF ( ip == -99999 )  THEN
 !--              wall category not found
-                 WRITE (message_string, "(A,I7,A,3I5)") 'wall category ', it,  &
-                                        ' not found  for i,j,k=', iw,jw,kw
-                 WRITE(9,*) message_string
+                 WRITE (9, "(A,I7,A,3I5)") 'wall category ', it,  &
+                                           ' not found  for i,j,k=', iw,jw,kw
+                 FLUSH(9)
+                 category = wall_category 
+                 DO k = 1, n_surface_types
+                    IF ( surface_type_codes(k) == category ) THEN
+                       ip = k
+                       EXIT
+                    ENDIF
+                 ENDDO
+                 IF ( ip == -99999 )  THEN
+!--                 default wall category not found
+                    WRITE (9, "(A,I5,A,3I5)") 'Default wall category', category, ' not found!'
+                    FLUSH(9)
+                    ip = 1
+                 ENDIF
               ENDIF
+
 !
 !--           Albedo
@@ -7192,4 +7400,8 @@
         ENDDO
 
+        
+        WRITE(9,*) 'Urban surfaces read'
+        FLUSH(9)
+        
         CALL location_message( '    types and parameters of urban surfaces read', .TRUE. )
    
@@ -7473,5 +7685,5 @@
            IF ( humidity )  surf_usm_h%vpt_surface(m) =                        &
                                                    surf_usm_h%pt_surface(m)
-                       
+
 !--        calculate true tendency
            stend = ( t_surf_h_p(m) - t_surf_h(m) - dt_3d * tsc(3) *           &
@@ -7508,5 +7720,5 @@
 !
 !--        calculate fluxes
-!--        rad_net_l is never used!           
+!--        rad_net_l is never used!
            surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) +                           &
                                      surf_usm_h%frac(ind_veg_wall,m) *                   &
@@ -7522,7 +7734,7 @@
            surf_usm_h%wghf_eb(m)   = lambda_surface *                         &
                                       ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) )
-           surf_usm_h%wghf_eb_green(m)  = lambda_surface_green *                        &
+           surf_usm_h%wghf_eb_green(m)  = lambda_surface_green *                         &
                                           ( t_surf_green_h_p(m) - t_green_h(nzb_wall,m) )
-           surf_usm_h%wghf_eb_window(m) = lambda_surface_window *                       &
+           surf_usm_h%wghf_eb_window(m) = lambda_surface_window *                        &
                                            ( t_surf_window_h_p(m) - t_window_h(nzb_wall,m) )
 
@@ -7539,5 +7751,5 @@
 !--        diffusion_s, surface_layer_fluxes,...
            surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / c_p
-
+      
        ENDDO
 !
@@ -7562,5 +7774,5 @@
 #if ! defined( __nopointer )          
 !
-!--          calculate rho * c_p coefficient at surface layer
+!--          calculate rho * c_p coefficient at wall layer
              rho_cp  = c_p * hyp(k) / ( r_d * surf_usm_v(l)%pt1(m) * exner(k) )
 #endif
@@ -7589,8 +7801,8 @@
 !--          obtained by simple linear interpolation. ( An alternative would
 !--          be an logarithmic interpolation. )
-!--          A roughness lenght of 0.001 is assumed for concrete (the inverse,
-!--          1000 is used in the nominator for scaling)
-             surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) * 1000.0_wp &
-                        * ( 11.8_wp + 4.2_wp *                                 &
+!--          Parameter roughness_concrete (default value = 0.001) is used
+!--          to calculation of roughness relative to concrete
+             surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) /           &
+                        roughness_concrete * ( 11.8_wp + 4.2_wp *              &
                         SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
                                    ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
@@ -7606,6 +7818,4 @@
              f_shf_window  = rho_cp / surf_usm_v(l)%r_a(m)
              f_shf_green   = rho_cp / surf_usm_v(l)%r_a(m)
-
-
 
 !--          add LW up so that it can be removed in prognostic equation
@@ -7657,6 +7867,4 @@
                            ( surf_usm_v(l)%c_surface_green(m) + coef_green_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) *       &
@@ -7677,5 +7885,5 @@
               IF ( humidity )  surf_usm_v(l)%vpt_surface(m) =                  &
                                                     surf_usm_v(l)%pt_surface(m)
-            
+
 !--           calculate true tendency
               stend = ( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) - dt_3d * tsc(3) *&
@@ -7766,20 +7974,21 @@
            dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
            dhour = INT(dtime/3600.0_wp)
-           DO m = 1, naheatlayers
-!--           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 ( k > get_topography_top_index_ji( j, i, 's' )  .AND. &
-                   k <= naheatlayers )  THEN
-!--              increase of pt in box i,j,k in time dt_3d 
-!--              given to anthropogenic heat aheat*acoef (W*m-2)
-!--              linear interpolation of coeficient
-                 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + &
-                         (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1)
-                 IF ( aheat(k,j,i) > 0.0_wp )  THEN
-                    pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exner(k)*rho_cp*dzu(k))
-                 ENDIF
-              ENDIF
+           DO i = nxl, nxr
+              DO j = nys, nyn
+                 DO k = nzub, min(nzut,naheatlayers)
+                    IF ( k > get_topography_top_index_ji( j, i, 's' ) ) THEN
+!--                    increase of pt in box i,j,k in time dt_3d 
+!--                    given to anthropogenic heat aheat*acoef (W*m-2)
+!--                    linear interpolation of coeficient
+                       acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + &
+                               (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1)
+                       IF ( aheat(k,j,i) > 0.0_wp )  THEN
+!--                       calculate rho * c_p coefficient at layer k
+                          rho_cp  = c_p * hyp(k) / ( r_d * pt(k+1,j,i) * exner(k) )
+                          pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exner(k)*rho_cp*dz(1))
+                       ENDIF
+                    ENDIF
+                 ENDDO
+              ENDDO
            ENDDO
 
