Changeset 2024 for palm/trunk

Timestamp:

Oct 12, 2016 4:42:37 PM (8 years ago)

Author:

kanani

Message:

changes related to urban surface model and output of ssws

Location:

palm/trunk/SOURCE

Files:

• check_parameters.f90 (modified) (5 diffs)
• plant_canopy_model_mod.f90 (modified) (2 diffs)
• sum_up_3d_data.f90 (modified) (2 diffs)
• urban_surface_mod.f90 (modified) (68 diffs)

palm/trunk/SOURCE/check_parameters.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added missing CASE, error message and unit for ssws*,
+! increased number of possible output quantities to 500.
 ! 
 ! Former revisions:
@@ -2834 +2835 @@
     IF ( data_output_user(1) /= ' ' )  THEN
        i = 1
-       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
+       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
           i = i + 1
        ENDDO
        j = 1
-       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
-          IF ( i > 100 )  THEN
+       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
+          IF ( i > 500 )  THEN
              message_string = 'number of output quantitities given by data' // &
-                '_output and data_output_user exceeds the limit of 100'
+                '_output and data_output_user exceeds the limit of 500'
              CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
           ENDIF
@@ -3018 +3019 @@
              CONTINUE
 
-          CASE ( 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 't*', &
+          CASE ( 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 'ssws*', 't*', &
                  'u*', 'z0*', 'z0h*', 'z0q*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
@@ -3054 +3055 @@
                 CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
              ENDIF
+             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res passive_scalar = .TRUE.'
+                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
+             ENDIF             
 
              IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
@@ -3061 +3067 @@
              IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
              IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
+             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'      
              IF ( TRIM( var ) == 't*'     )  unit = 'K'
              IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'

palm/trunk/SOURCE/plant_canopy_model_mod.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added missing lad_s initialization
 ! 
 ! Former revisions:
@@ -707 +707 @@
         REAL(wp), DIMENSION(:), ALLOCATABLE     :: col
 
+        lad_s = 0.0_wp
         OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
                 action='READ', status='OLD', form='FORMATTED', err=515)

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added missing CASE for ssws*
 ! 
 ! Former revisions:
@@ -343 +343 @@
                 ENDIF
                 shf_av = 0.0_wp
+                
+             CASE ( 'ssws*' )
+                IF ( .NOT. ALLOCATED( ssws_av ) )  THEN
+                   ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                ssws_av = 0.0_wp                
 
              CASE ( 't*' )

palm/trunk/SOURCE/urban_surface_mod.f90

@@ -21 +21 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
+! optimization of MPI-RMA operations,
+! declaration of pcbl as integer,
+! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
+! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
+! use of new control parameter varnamelength,
+! added output variables usm_rad_ressw, usm_rad_reslw,
+! minor formatting changes,
+! minor optimizations.
 ! 
 ! Former revisions:
@@ -43 +51 @@
 ! ------------
 ! 2016/6/9 - Initial version of the USM (Urban Surface Model)
-!            authors: Jaroslav Resler, Pavel Krc (CTU in Prague, ICS AS in Prague)
+!            authors: Jaroslav Resler, Pavel Krc
+!                     (Czech Technical University in Prague and Institute of
+!                      Computer Science of the Czech Academy of Sciences, Prague)
 !            with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
 !            partly inspired by PALM LSM (B. Maronga)
@@ -73 +83 @@
 !>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
 !>
+!> @todo Check optimizations for RMA operations
+!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
+!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
+!>       factor 3 between min and max time
 !------------------------------------------------------------------------------!
  MODULE urban_surface_mod
@@ -92 +106 @@
                time_since_reference_point, surface_pressure,                   &
                g, pt_surface, large_scale_forcing, lsf_surf,                   &
-               time_do3d, dt_do3d, average_count_3d, urban_surface
+               time_do3d, dt_do3d, average_count_3d, varnamelength,            &
+               urban_surface
 
     USE cpulog,                                                                &
@@ -161 +176 @@
     INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers in urban surface layer above top of buildings 
     INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
+    INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
     INTEGER(iwp), PARAMETER                        ::  ndcsf = 2                          !< number of dimensions of real values in CSF
-    INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF
+    INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
+    INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
     INTEGER(iwp), PARAMETER                        ::  id = 1                             !< position of d-index in surfl and surf
     INTEGER(iwp), PARAMETER                        ::  iz = 2                             !< position of k-index in surfl and surf
@@ -205 +222 @@
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf
                                                                         !< respective block for particular processor is surfstart[iproc]+1 : surfstart[iproc+1]
-    INTEGER(iwp)                                   ::  nsvfl            !< number of svf (excluding csf) for local processor
-    INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor 
+    INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
+    INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
                                                                         !< needed only during calc_svf but must be here because it is
-                                                                        !< shared between subroutines usm_calc_svf and usm_raytrace                                                                         
-    INTEGER(iwp)                                   ::  nsvfcsfl         !< sum of svf+csf for local processor
-
+                                                                        !< shared between subroutines usm_calc_svf and usm_raytrace
 
 !-- type for calculation of svf
@@ -240 +255 @@
     INTEGER(iwp)                                   ::  msvf, mcsf       !< mod for swapping the growing array
     INTEGER(iwp), PARAMETER                        ::  gasize = 10000   !< initial size of growing arrays
+!-- temporary arrays for calculation of csf in raytracing
+    INTEGER(iwp)                                   ::  maxboxesg        !< max number of boxes ray can cross in the domain
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  boxes            !< coordinates of gridboxes being crossed by ray
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  crlens           !< array of crossing lengths of ray for particular grid boxes
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored 
+    INTEGER(kind=MPI_ADDRESS_KIND), &
+                  DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray 
 
 !-- arrays storing the values of USM
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of source and target surface for svf[isvf]
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors
-                                                                        !< for individual local surfaces and plant canopy sinks
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl          !< array of lw radiation for local surface after i-th reflection
@@ -276 +298 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
+    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    
     
 !-- block variables needed for calculation of the plant canopy model inside the urban surface model
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  pcbl             !< z,y,x coordinates of i-th local plant canopy box pcbl[:,i] = [z, y, x]
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< index of local pcb[z,y,x]
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  csfsurf          !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  csf              !< array of plant canopy sink fators + direct irradiation factors (transparency)
+                                                                        !< for local surfaces
+    INTEGER(wp), DIMENSION(:,:), ALLOCATABLE       ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< index of local pcb[k,j,i]
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box
@@ -486 +513 @@
         INTEGER(iwp)                            :: i, j, k, d, l, ir, jr, ids
         INTEGER(iwp)                            :: nzubl, nzutl, isurf, ipcgb
-        INTEGER                                 :: procid
+        INTEGER(iwp)                            :: procid
 
         
@@ -828 +855 @@
   
         INTEGER(iwp)                                       :: i, j, k, l, ids, iwl,istat
-        CHARACTER (len=20)                                 :: var, surfid
+        CHARACTER (len=varnamelength)                      :: var, surfid
         INTEGER(iwp), PARAMETER                            :: nd = 5
         CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER     :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
@@ -861 +888 @@
            SELECT CASE ( TRIM( var ) )
                 
-                CASE ( 'usm_radnet' )
+                CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
                     IF ( .NOT.  ALLOCATED(rad_net_av) )  THEN
@@ -931 +958 @@
                     ENDIF
 
+                CASE ( 'usm_rad_ressw' )
+!--                 array of residua of sw radiation absorbed in surface after last reflection
+                    IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
+                        ALLOCATE( surfins_av(startenergy:endenergy) )
+                        surfins_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_reslw' )
+!--                 array of residua of lw radiation absorbed in surface after last reflection
+                    IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
+                        ALLOCATE( surfinl_av(startenergy:endenergy) )
+                        surfinl_av = 0.0_wp
+                    ENDIF
+                                    
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
@@ -959 +1000 @@
                         t_surf_av = 0.0_wp
                     ENDIF
-                    
+
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for iwl layer of walls and land
@@ -976 +1017 @@
            SELECT CASE ( TRIM( var ) )
                 
-                CASE ( 'usm_radnet' )
+                CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
                     DO l = startenergy, endenergy
@@ -1058 +1099 @@
                     ENDDO
                     
+                CASE ( 'usm_rad_ressw' )
+!--                 array of residua of sw radiation absorbed in surface after last reflection
+                    DO l = startenergy, endenergy
+                        IF ( surfl(id,l) == ids )  THEN
+                            surfins_av(l) = surfins_av(l) + surfins(l)
+                        ENDIF
+                    ENDDO
+                                    
+                CASE ( 'usm_rad_reslw' )
+!--                 array of residua of lw radiation absorbed in surface after last reflection
+                    DO l = startenergy, endenergy
+                        IF ( surfl(id,l) == ids )  THEN
+                            surfinl_av(l) = surfinl_av(l) + surfinl(l)
+                        ENDIF
+                    ENDDO
+                    
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
@@ -1107 +1164 @@
            SELECT CASE ( TRIM( var ) )
                 
-                CASE ( 'usm_radnet' )
+                CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
                     DO l = startenergy, endenergy
@@ -1187 +1244 @@
                     ENDDO
 
+                CASE ( 'usm_rad_ressw' )
+!--                 array of residua of sw radiation absorbed in surface after last reflection
+                    DO l = startenergy, endenergy
+                        IF ( surfl(id,l) == ids )  THEN
+                            surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
+                        ENDIF
+                    ENDDO
+                                    
+                CASE ( 'usm_rad_reslw' )
+!--                 array of residua of lw radiation absorbed in surface after last reflection
+                    DO l = startenergy, endenergy
+                        IF ( surfl(id,l) == ids )  THEN
+                            surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
+                        ENDIF
+                    ENDDO
+                    
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
@@ -1218 +1291 @@
                         ENDIF
                     ENDDO
-                    
+
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for  iwl layer of walls and land
@@ -1226 +1299 @@
                         ENDIF
                     ENDDO
-                
+
            END SELECT
 
@@ -1391 +1464 @@
         REAL(wp), DIMENSION(3)                      :: uv
         LOGICAL                                     :: visible
-        REAL(wp)                                    :: sz, sy, sx, tz, ty, tx, transparency, rirrf, sqdist, svfsum
-        !REAL(wp)                                    :: rsvf
+        REAL(wp), DIMENSION(3)                      :: sa, ta          !< real coordinates z,y,x of source and target
+        REAL(wp)                                    :: transparency, rirrf, sqdist, svfsum
         INTEGER(iwp)                                :: isurflt, isurfs, isurflt_prev
         INTEGER(iwp)                                :: itx, ity, itz
@@ -1439 +1512 @@
         IF ( plant_canopy )  THEN
             ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) )
+            maxboxesg = nx + ny + nzu + 1
+!--         temporary arrays storing values for csf calculation during raytracing
+            ALLOCATE( boxes(3, maxboxesg) )
+            ALLOCATE( crlens(maxboxesg) )
+
 #if defined( __parallel )
             ALLOCATE( planthl(nys:nyn,nxl:nxr) )
@@ -1445 +1523 @@
             CALL MPI_AllGather( planthl, nnx*nny, MPI_INTEGER, &
                                 plantt, nnx*nny, MPI_INTEGER, comm2d, ierr )
-            DEALLOCATE(planthl)
+            DEALLOCATE( planthl )
+            
+!--         temporary arrays storing values for csf calculation during raytracing
+            ALLOCATE( lad_ip(maxboxesg) )
+            ALLOCATE( lad_disp(maxboxesg) )
 
             IF ( usm_lad_rma )  THEN
+                ALLOCATE( lad_s_ray(maxboxesg) )
+                
+                ! set conditions for RMA communication
                 CALL MPI_Info_create(minfo, ierr)
                 CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr)
@@ -1500 +1585 @@
                 IF ( i < nxl  .OR.  i > nxr &
                      .OR.  j < nys  .OR.  j > nyn ) CYCLE
-                tx = REAL(i)
-                ty = REAL(j)
-                tz = REAL(k)
+                ta = (/ REAL(k), REAL(j), REAL(i) /)
 
                 DO isurfs = 1, nsurf
@@ -1512 +1595 @@
                       
                     sd = surf(id, isurfs)
-                    sz = REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd)
-                    sy = REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd)
-                    sx = REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)
+                    sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &
+                            REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &
+                            REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)
 
 !--                 unit vector source -> target
-                    uv = (/ (tz-sz)*dz, (ty-sy)*dy, (tx-sx)*dx /)
+                    uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
                     sqdist = SUM(uv(:)**2)
                     uv = uv / SQRT(sqdist)
@@ -1527 +1610 @@
 
 !--                 raytrace while not creating any canopy sink factors
-                    CALL usm_raytrace((/sz,sy,sx/), (/tz,ty,tx/), isurfs, rirrf, 1._wp, .FALSE., &
+                    CALL usm_raytrace(sa, ta, isurfs, rirrf, 1._wp, .FALSE., &
                             visible, transparency, win_lad)
                     IF ( .NOT.  visible ) CYCLE
@@ -1558 +1641 @@
             td = surfl(id, isurflt)
             IF ( td >= isky  .AND.  .NOT.  plant_canopy ) CYCLE
-            tz = REAL(surfl(iz, isurflt), wp) - 0.5_wp * kdir(td)
-            ty = REAL(surfl(iy, isurflt), wp) - 0.5_wp * jdir(td)
-            tx = REAL(surfl(ix, isurflt), wp) - 0.5_wp * idir(td)
+            ta = (/ REAL(surfl(iz, isurflt), wp) - 0.5_wp * kdir(td),  &
+                      REAL(surfl(iy, isurflt), wp) - 0.5_wp * jdir(td),  &
+                      REAL(surfl(ix, isurflt), wp) - 0.5_wp * idir(td)  /)
             DO isurfs = 1, nsurf
                 IF ( .NOT.  usm_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
@@ -1570 +1653 @@
                   
                 sd = surf(id, isurfs)
-                sz = REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd)
-                sy = REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd)
-                sx = REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)
+                sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
+                        REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
+                        REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
 
 !--             unit vector source -> target
-                uv = (/ (tz-sz)*dz, (ty-sy)*dy, (tx-sx)*dx /)
+                uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
                 sqdist = SUM(uv(:)**2)
                 uv = uv / SQRT(sqdist)
@@ -1586 +1669 @@
 
 !--             raytrace + process plant canopy sinks within
-                CALL usm_raytrace((/sz,sy,sx/), (/tz,ty,tx/), isurfs, rirrf, facearea(td), .TRUE., &
+                CALL usm_raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
                         visible, transparency, win_lad)
                 
@@ -1623 +1706 @@
 
         CALL location_message( '    waiting for completion of SVF and CSF calculation in all processes', .TRUE. )
-
+!--     deallocate temporary global arrays
+        DEALLOCATE(nzterr)
+        
+        IF ( plant_canopy )  THEN
+!--         finalize mpi_rma communication and deallocate temporary arrays
 #if defined( __parallel )
-        IF ( plant_canopy )  THEN
             IF ( usm_lad_rma )  THEN
                 CALL MPI_Win_flush_all(win_lad, ierr)
@@ -1632 +1718 @@
 !--             free MPI window
                 CALL MPI_Win_free(win_lad, ierr)
+                
+!--             deallocate temporary arrays storing values for csf calculation during raytracing
+                DEALLOCATE( lad_s_ray )
+!--             usm_lad is the pointer to lad_s_rma in case of usm_lad_rma
+!--             and must not be deallocated here
             ELSE
                 DEALLOCATE(usm_lad)
                 DEALLOCATE(usm_lad_g)
             ENDIF
+#else
+            DEALLOCATE(usm_lad)
+#endif
+            DEALLOCATE( boxes )
+            DEALLOCATE( crlens )
+            DEALLOCATE( plantt )
         ENDIF
-#else
-        DEALLOCATE(usm_lad)
-#endif
-
-!--     deallocate temporary global arrays
-        IF ( ALLOCATED(nzterr) ) DEALLOCATE(nzterr)
-        IF ( ALLOCATED(plantt) ) DEALLOCATE(plantt)
-        
+
+        CALL location_message( '    calculation of the complete SVF array', .TRUE. )
+
 !--     sort svf ( a version of quicksort )
         CALL quicksort_svf(asvf,1,nsvfl)
-        
+
+        ALLOCATE( svf(ndsvf,nsvfl) )
+        ALLOCATE( svfsurf(idsvf,nsvfl) )
+
+        !< load svf from the structure array to plain arrays
+        isurflt_prev = -1
+        ksvf = 1
+        svfsum = 0._wp
+        DO isvf = 1, nsvfl
+!--         normalize svf per target face
+            IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
+                IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
+!--                 TODO detect and log when normalization differs too much from 1
+                    svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum
+                ENDIF
+                isurflt_prev = asvf(ksvf)%isurflt
+                isvf_surflt = isvf
+                svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
+            ELSE
+                svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
+            ENDIF
+
+            svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
+            svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
+
+!--         next element
+            ksvf = ksvf + 1
+        ENDDO
+
+        IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
+!--         TODO detect and log when normalization differs too much from 1
+            svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum
+        ENDIF
+
+!--     deallocate temporary asvf array
+!--     DEALLOCATE(asvf) - ifort has a problem with deallocation of allocatable target
+!--     via pointing pointer - we need to test original targets
+        IF ( ALLOCATED(asvf1) )  THEN
+            DEALLOCATE(asvf1)
+        ENDIF
+        IF ( ALLOCATED(asvf2) )  THEN
+            DEALLOCATE(asvf2)
+        ENDIF
+
         npcsfl = 0
         IF ( plant_canopy )  THEN
+
+            CALL location_message( '    calculation of the complete CSF array', .TRUE. )
+
 !--         sort and merge csf for the last time, keeping the array size to minimum
             CALL usm_merge_and_grow_csf(-1)
@@ -1720 +1858 @@
 !--         scatter and gather the number of elements to and from all processor
 !--         and calculate displacements
-            CALL mpi_alltoall(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
+            CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
             
             npcsfl = SUM(ipcsflt)
@@ -1788 +1926 @@
                 npcsfl = kcsf
             ENDIF
-            
-        ENDIF  !< plant_canopy
-        
-        CALL location_message( '    calculation of the complete SVF array', .TRUE. )
-        nsvfcsfl = nsvfl + npcsfl
-        
-        ALLOCATE( svf(ndsvf,nsvfcsfl) )
-        ALLOCATE( svfsurf(ndsvf,nsvfcsfl) )
-
-        !< load svf from the structure array to plain arrays
-        isurflt_prev = -1
-        ksvf = 1
-        svfsum = 0._wp
-        DO isvf = 1, nsvfl
-!--         normalize svf per target face
-            IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
-                IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
-!--                 TODO detect and log when normalization differs too much from 1
-                    svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum
-                ENDIF
-                isurflt_prev = asvf(ksvf)%isurflt
-                isvf_surflt = isvf
-                svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
-            ELSE
-                svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
-            ENDIF
-            
-            svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
-            svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
-                    
-!--         next element
-            ksvf = ksvf + 1    
-        ENDDO
-        
-        IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
-!--         TODO detect and log when normalization differs too much from 1
-            svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum
-        ENDIF
-        
-!--     deallocate temporary asvf array
-!--     DEALLOCATE(asvf) - ifort has a problem with deallocation of allocatable target
-!--     via pointing pointer - we need to test original targets
-        IF ( ALLOCATED(asvf1) )  THEN
-            DEALLOCATE(asvf1)
-        ENDIF
-        IF ( ALLOCATED(asvf2) )  THEN
-            DEALLOCATE(asvf2)
-        ENDIF
-
-        IF ( plant_canopy )  THEN
-            CALL location_message( '    calculation of the complete CSF part of the array', .TRUE. )
-            IF ( npcsfl > 0 )  THEN
-                DO isvf = 1, npcsfl
-                    svf(:,nsvfl+isvf) = pcsflt(:,isvf)
-                    svfsurf(1,nsvfl+isvf) =  gridpcbl(kpcsflt(1,isvf),kpcsflt(2,isvf),kpcsflt(3,isvf))
-                    svfsurf(2,nsvfl+isvf) =  kpcsflt(4,isvf)
+
+            ncsfl = npcsfl
+            IF ( ncsfl > 0 )  THEN
+                ALLOCATE( csf(ndcsf,ncsfl) )
+                ALLOCATE( csfsurf(idcsf,ncsfl) )
+                DO icsf = 1, ncsfl
+                    csf(:,icsf) = pcsflt(:,icsf)
+                    csfsurf(1,icsf) =  gridpcbl(kpcsflt(1,icsf),kpcsflt(2,icsf),kpcsflt(3,icsf))
+                    csfsurf(2,icsf) =  kpcsflt(4,icsf)
                 ENDDO
             ENDIF
@@ -1877 +1968 @@
         CHARACTER (len=*),INTENT(OUT)   ::  unit     !:
         
-        CHARACTER (len=20)              :: var
+        CHARACTER (len=varnamelength)   :: var
 
         var = TRIM(variable)
-        IF ( var(1:11)=='usm_radnet_'  .OR.  var(1:13)=='usm_rad_insw_'  .OR.           &
-             var(1:13)=='usm_rad_inlw_'  .OR.  var(1:16)=='usm_rad_inswdir_'  .OR.      &
-             var(1:16)=='usm_rad_inswdif_'  .OR.  var(1:16)=='usm_rad_inswref_'  .OR.   &
-             var(1:16)=='usm_rad_inlwdif_'  .OR.  var(1:16)=='usm_rad_inlwref_'  .OR.   &
-             var(1:14)=='usm_rad_outsw_'  .OR.  var(1:14)=='usm_rad_outlw_'  .OR.       &
-             var(1:11)=='usm_rad_hf_'  .OR.                                             &
-             var(1:9) =='usm_wshf_'  .OR.  var(1:9)=='usm_wghf_' )  THEN
+        IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.        &
+             var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.    &
+             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR. &
+             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR. &
+             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:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' )  THEN
             unit = 'W/m2'
-        ELSE IF ( var(1:10) =='usm_t_surf'  .OR.  var(1:10) =='usm_t_wall' )  THEN
+        ELSE IF ( var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall' )  THEN
             unit = 'K'
-        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.           &
-                  var(1:11) =='usm_surfalb'  .OR.  var(1:12) =='usm_surfemis')  THEN
+        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.           &
+                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis')  THEN
             unit = '1'
-        ELSE IF ( plant_canopy  .AND.  var(1:7) =='usm_lad' )  THEN
+        ELSE IF ( plant_canopy  .AND.  var(1:7) == 'usm_lad' )  THEN
             unit = 'm2/m3'
-        ELSE IF ( plant_canopy  .AND.  var(1:14) == 'usm_canopy_khf' )  THEN
+        ELSE IF ( plant_canopy  .AND.  var(1:13) == 'usm_canopy_hr' )  THEN
             unit = 'K/s'
         ELSE
@@ -1972 +2064 @@
         REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)     ::  temp_pf    !< temp array for urban surface output procedure
         
-        CHARACTER (len=20)                                     :: var, surfid
+        CHARACTER (len=varnamelength)                          :: var, surfid
         INTEGER(iwp), PARAMETER                                :: nd = 5
         CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER         :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
@@ -2085 +2177 @@
               ENDDO
 
-          CASE ( 'usm_radnet' )
+          CASE ( 'usm_rad_net' )
 !--           array of complete radiation balance
               DO isurf = dirstart(ids), dirend(ids)
@@ -2202 +2294 @@
                    ELSE
                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
+                   ENDIF
+                 ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_ressw' )
+!--           average of array of residua of sw radiation absorbed in surface after last 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)) = surfins(isurf)
+                   ELSE
+                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
+                   ENDIF
+                 ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_reslw' )
+!--           average of array of residua of lw radiation absorbed in surface after last 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)) = surfinl(isurf)
+                   ELSE
+                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
                    ENDIF
                  ENDIF
@@ -2277 +2393 @@
               ENDDO
               
-          CASE ( 'usm_canopy_khf' )
-!--           canopy kinematic heat flux
+          CASE ( 'usm_canopy_hr' )
+!--           canopy heating rate
               DO i = nxl, nxr
                  DO j = nys, nyn
@@ -2322 +2438 @@
         CHARACTER (len=*), INTENT(OUT) ::  grid_z      !< 
 
-        CHARACTER (len=20)            :: var
+        CHARACTER (len=varnamelength)  :: var
 
         var = TRIM(variable)
-        IF ( var(1:11)=='usm_radnet_'  .OR.  var(1:13) =='usm_rad_insw_'  .OR.            &
-             var(1:13) =='usm_rad_inlw_'  .OR.  var(1:16) =='usm_rad_inswdir_'  .OR.      &
-             var(1:16) =='usm_rad_inswdif_'  .OR.  var(1:16) =='usm_rad_inswref_'  .OR.   &
-             var(1:16) =='usm_rad_inlwdif_'  .OR.  var(1:16) =='usm_rad_inlwref_'  .OR.   &
-             var(1:14) =='usm_rad_outsw_'  .OR.  var(1:14) =='usm_rad_outlw_'  .OR.       &
-             var(1:11) =='usm_rad_hf_'  .OR.                                              &
-             var(1:9) == 'usm_wshf_'  .OR.  var(1:9)== 'usm_wghf_'  .OR.                  &
-             var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall'  .OR.             &
-             var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                   & 
-             var(1:7) =='usm_dif'  .OR.  var(1:11) =='usm_surfcat'  .OR.                  &
-             var(1:11) =='usm_surfalb'  .OR.  var(1:12) =='usm_surfemis'  .OR.            &
-             var(1:7) == 'usm_lad'  .OR.  var(1:14) == 'usm_canopy_khf' )  THEN
+        IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.          &
+             var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.      &
+             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR.   &
+             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR.   &
+             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:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
+             var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall'  .OR.               &
+             var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                     & 
+             var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
+             var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.            &
+             var(1:7) == 'usm_lad'  .OR.  var(1:13) == 'usm_canopy_hr' )  THEN
 
             found = .TRUE.
@@ -2531 +2648 @@
         
         IF ( read_svf_on_init )  THEN
-!--         read svf and svfsurf data from file
+!--         read svf, csf, svfsurf and csfsurf data from file
             CALL location_message( '    Start reading SVF from file', .TRUE. )
             CALL usm_read_svf_from_file()
@@ -2545 +2662 @@
 
         IF ( write_svf_on_init )  THEN
-!--         write svf and svfsurf data to file
+!--         write svf, csf svfsurf and csfsurf data to file
+            CALL location_message( '    Store SVF and CSF to file', .TRUE. )
             CALL usm_write_svf_to_file()
         ENDIF
@@ -2686 +2804 @@
 ! Description:
 ! ------------
+!> Parin for &usm_par for urban surface model
+!------------------------------------------------------------------------------!
+    SUBROUTINE usm_parin
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=80) ::  line  !< string containing current line of file PARIN
+
+       NAMELIST /urban_surface_par/                                            &
+                           land_category,                                      &
+                           mrt_factors,                                        &
+                           nrefsteps,                                          &
+                           pedestrant_category,                                &
+                           ra_horiz_coef,                                      &
+                           read_svf_on_init,                                   &
+                           roof_category,                                      &
+                           split_diffusion_radiation,                          &
+                           urban_surface,                                      &
+                           usm_anthropogenic_heat,                             &
+                           usm_energy_balance_land,                            &
+                           usm_energy_balance_wall,                            &
+                           usm_material_model,                                 &
+                           usm_lad_rma,                                        &
+                           wall_category,                                      &
+                           write_svf_on_init
+
+       line = ' '
+
+!
+!--    Try to find urban surface model package
+       REWIND ( 11 )
+       line = ' '
+       DO   WHILE ( INDEX( line, '&urban_surface_par' ) == 0 )
+          READ ( 11, '(A)', END=10 )  line
+       ENDDO
+       BACKSPACE ( 11 )
+
+!
+!--    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
+
+    END SUBROUTINE usm_parin
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> This subroutine calculates interaction of the solar radiation
 !> with urban surface and updates surface, roofs and walls heatfluxes.
@@ -2695 +2867 @@
        
         INTEGER(iwp)               :: i, j, k, kk, is, js, d, ku, refstep
-        INTEGER(iwp)               :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, ipcgb
+        INTEGER(iwp)               :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb
         INTEGER(iwp), DIMENSION(4) :: bdycross
         REAL(wp), DIMENSION(3,3)   :: mrot            !< grid rotation matrix (xyz)
@@ -2836 +3008 @@
 !--         pcsf first pass
             isurf1 = -1  !< previous processed pcgb
-            DO isvf = nsvfl+1, nsvfcsfl
-                ipcgb = svfsurf(1, isvf)
+            DO icsf = 1, ncsfl
+                ipcgb = csfsurf(1, icsf)
                 i = pcbl(ix,ipcgb)
                 j = pcbl(iy,ipcgb)
                 k = pcbl(iz,ipcgb)
-                isurfsrc = svfsurf(2, isvf)
+                isurfsrc = csfsurf(2, icsf)
 
                 IF ( zenith(0) > 0  .AND.  ipcgb /= isurf1 )  THEN
@@ -2918 +3090 @@
 
 !--         radiation absorbed by plant canopy
-            DO isvf = nsvfl+1, nsvfcsfl
-                ipcgb = svfsurf(1, isvf)
-                isurfsrc = svfsurf(2, isvf)
+            DO icsf = 1, ncsfl
+                ipcgb = csfsurf(1, icsf)
+                isurfsrc = csfsurf(2, icsf)
 
                 IF ( surf(id, isurfsrc) < isky )  THEN
-                    pcbinsw(ipcgb) = pcbinsw(ipcgb) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc)
-!--                 pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc)
+                    pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
+!--                 pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc)
                 ENDIF
             ENDDO
@@ -2998 +3170 @@
         REAL(wp), INTENT(out)                  :: transparency !< along whole path
         INTEGER(iwp), INTENT(in)               :: win_lad
-        INTEGER(iwp)                           :: k, d
+        INTEGER(iwp)                           :: i, j, k, d
         INTEGER(iwp)                           :: seldim       !< dimension to be incremented
         INTEGER(iwp)                           :: ncsb         !< no of written plant canopy sinkboxes
@@ -3017 +3189 @@
         INTEGER(iwp)                           :: px, py       !< number of processors in x and y dir before 
                                                                !< the processor in the question
-        INTEGER(iwp)                           :: ig, ip
-        REAL(wp)                               :: lad_s_target
-        INTEGER(kind=MPI_ADDRESS_KIND)         :: lad_disp
-        REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp
+        INTEGER(iwp)                           :: ip           !< number of processor where gridbox reside
+        INTEGER(iwp)                           :: ig           !< 1D index of gridbox in global 2D array
+        REAL(wp)                               :: lad_s_target !< recieved lad_s of particular grid box
+        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
@@ -3086 +3258 @@
                 IF ( box(1) <= nzterr(ig) )  THEN
                     visible = .FALSE.
-                    IF ( ncsb > 0 )  THEN
-!--                     rewind written plant canopy sink factors - they are invalid
-                        ncsfl = ncsfl - ncsb
-                    ENDIF
                     RETURN
                 ENDIF
@@ -3095 +3263 @@
                 IF ( plant_canopy )  THEN
                     IF ( box(1) <= plantt(ig) )  THEN
+                        ncsb = ncsb + 1
+                        boxes(:,ncsb) = box
+                        crlens(ncsb) = crlen
 #if defined( __parallel )
-                        lad_disp = (box(3)-px*nnx)*(nny*nzu) + (box(2)-py*nny)*nzu + box(1)-nzub
-                        IF ( usm_lad_rma )  THEN
-!--                         Read LAD using MPI RMA
-                            CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
-                            CALL MPI_Get(lad_s_target, 1, MPI_REAL, ip, lad_disp, 1, MPI_REAL, &
-                                win_lad, ierr)
-                            IF ( ierr /= 0 )  THEN
-                                WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
-                                CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
-                            ENDIF
-                            
-                            CALL MPI_Win_flush_local(ip, win_lad, ierr)
-                            IF ( ierr /= 0 )  THEN
-                                WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local'
-                                CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
-                            ENDIF
-                            CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
-                        ELSE
-                            lad_s_target = usm_lad_g(ip*nnx*nny*nzu + lad_disp)
-                        ENDIF
-#else
-                        lad_s_target = usm_lad(box(1),box(2),box(3))
+                        lad_ip(ncsb) = ip
+                        lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nzu) + (box(2)-py*nny)*nzu + box(1)-nzub
 #endif
-                        cursink = 1._wp - exp(-ext_coef * lad_s_target &
-                            * crlen*realdist)
-
-                        IF ( create_csf )  THEN
-!--                         write svf values into the array
-                            ncsb = ncsb + 1
-                            ncsfl = ncsfl + 1
-                            acsf(ncsfl)%ip = ip
-                            acsf(ncsfl)%itx = box(3)
-                            acsf(ncsfl)%ity = box(2)
-                            acsf(ncsfl)%itz = box(1)
-                            acsf(ncsfl)%isurfs = isrc
-                            acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !-- we postpone multiplication by transparency
-                            acsf(ncsfl)%rtransp = REAL(transparency, wp)
-                        ENDIF  !< create_csf
-
-                        transparency = transparency * (1._wp - cursink)
-                        
                     ENDIF
                 ENDIF
@@ -3147 +3280 @@
         ENDDO
         
+        IF ( plant_canopy )  THEN
+#if defined( __parallel )
+            IF ( usm_lad_rma )  THEN
+!--             send requests for lad_s to appropriate processor
+                CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
+                DO i = 1, ncsb
+                    CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), &
+                                 1, MPI_REAL, win_lad, ierr)
+                    IF ( ierr /= 0 )  THEN
+                        WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
+                        CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
+                    ENDIF
+                ENDDO
+                
+!--             wait for all pending local requests complete
+                CALL MPI_Win_flush_local_all(win_lad, ierr)
+                IF ( ierr /= 0 )  THEN
+                    WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
+                    CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
+                ENDIF
+                CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
+                
+            ENDIF
+#endif
+
+!--         calculate csf and transparency
+            DO i = 1, ncsb
+#if defined( __parallel )
+                IF ( usm_lad_rma )  THEN
+                    lad_s_target = lad_s_ray(i)
+                ELSE
+                    lad_s_target = usm_lad_g(lad_ip(i)*nnx*nny*nzu + lad_disp(i))
+                ENDIF
+#else
+                lad_s_target = usm_lad(boxes(1,i),boxes(2,i),boxes(3,i))
+#endif
+                cursink = 1._wp - exp(-ext_coef * lad_s_target * crlens(i)*realdist)
+
+                IF ( create_csf )  THEN
+!--                 write svf values into the array
+                    ncsfl = ncsfl + 1
+                    acsf(ncsfl)%ip = lad_ip(i)
+                    acsf(ncsfl)%itx = boxes(3,i)
+                    acsf(ncsfl)%ity = boxes(2,i)
+                    acsf(ncsfl)%itz = boxes(1,i)
+                    acsf(ncsfl)%isurfs = isrc
+                    acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !-- we postpone multiplication by transparency
+                    acsf(ncsfl)%rtransp = REAL(transparency, wp)
+                ENDIF  !< create_csf
+
+                transparency = transparency * (1._wp - cursink)
+                
+            ENDDO
+        ENDIF
+        
         visible = .TRUE.
+        
     END SUBROUTINE usm_raytrace
     
@@ -3197 +3386 @@
             
 #if defined( __parallel ) && ! defined ( __check )
-            CALL mpi_barrier( comm2d, ierr )
+            CALL MPI_BARRIER( comm2d, ierr )
 #endif
         ENDDO
@@ -3231 +3420 @@
             
 #if defined( __parallel ) && ! defined ( __check )
-            CALL mpi_barrier( comm2d, ierr )
+            CALL MPI_BARRIER( comm2d, ierr )
 #endif
         ENDDO
@@ -3256 +3445 @@
        CHARACTER (LEN=30) ::  variable_chr  !< dummy variable to read string
        
-       INTEGER            ::  i             !< running index      
-       
+       INTEGER(iwp)       ::  i             !< running index
+
 
        DO  i = 0, io_blocks-1
@@ -3317 +3506 @@
 
         IMPLICIT NONE
-        INTEGER                      :: fsvf = 89
-        INTEGER                      :: i
+        INTEGER(iwp)                 :: fsvf = 89
+        INTEGER(iwp)                 :: i
         CHARACTER(usm_version_len)   :: usm_version_field
         CHARACTER(svf_code_len)      :: svf_code_field
@@ -3336 +3525 @@
                 ENDIF
                 
-!--             read nsvfcsfl, nsvfl
-                READ ( fsvf ) nsvfcsfl, nsvfl
-                IF ( nsvfcsfl <= 0 )  THEN
+!--             read nsvfl, ncsfl
+                READ ( fsvf ) nsvfl, ncsfl
+                IF ( nsvfl <= 0  .OR.  ncsfl < 0 )  THEN
                     WRITE( message_string, * ) 'Wrong number of SVF or CSF'
                     CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
                 ELSE
-                    WRITE(message_string,*) '    Number of SVF and CSF to read', nsvfcsfl, nsvfl
+                    WRITE(message_string,*) '    Number of SVF and CSF to read', nsvfl, ncsfl
                     CALL location_message( message_string, .TRUE. )
                 ENDIF
                 
-                ALLOCATE(svf(ndsvf,nsvfcsfl))
-                ALLOCATE(svfsurf(ndsvf,nsvfcsfl))
-                
+                ALLOCATE(svf(ndsvf,nsvfl))
+                ALLOCATE(svfsurf(idsvf,nsvfl))
                 READ(fsvf) svf
                 READ(fsvf) svfsurf
+                IF ( plant_canopy )  THEN
+                    ALLOCATE(csf(ndcsf,ncsfl))
+                    ALLOCATE(csfsurf(idcsf,ncsfl))
+                    READ(fsvf) csf
+                    READ(fsvf) csfsurf
+                ENDIF
                 READ ( fsvf ) svf_code_field
                 
@@ -3362 +3556 @@
             ENDIF
 #if defined( __parallel )
-            CALL mpi_barrier( comm2d, ierr )
+            CALL MPI_BARRIER( comm2d, ierr )
 #endif
         ENDDO
@@ -3379 +3573 @@
     
         CHARACTER(12)                                         :: wtn
-        INTEGER                                               :: wtc
+        INTEGER(iwp)                                          :: wtc
         REAL(wp), DIMENSION(n_surface_params)                 :: wtp
     
@@ -3490 +3684 @@
             ENDIF
 #if defined( __parallel ) && ! defined ( __check )
-            CALL mpi_barrier( comm2d, ierr )
+            CALL MPI_BARRIER( comm2d, ierr )
 #endif
         ENDDO
@@ -3912 +4106 @@
        THEN
 #if defined( __parallel )
-          IF ( collective_wait )  CALL mpi_barrier( comm2d, ierr )
+          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 )
@@ -3934 +4128 @@
        IMPLICIT NONE
 
-       INTEGER, INTENT(IN) :: mod_count
-       INTEGER :: i
+       INTEGER(iwp), INTENT(IN) :: mod_count
+       INTEGER(iwp)             :: i
       
 #if defined( __nopointer )
@@ -4027 +4221 @@
        IMPLICIT NONE
        
-       INTEGER ::  i
+       INTEGER(iwp)  ::  i
 
        DO  i = 0, io_blocks-1
@@ -4058 +4252 @@
 ! Description:
 ! ------------
-!> Subroutine stores svf and svfsurf data to files
+!> Subroutine stores svf, svfsurf, csf and csfsurf data to a file.
 !------------------------------------------------------------------------------!
     SUBROUTINE usm_write_svf_to_file
-    
+
         IMPLICIT NONE
-        INTEGER             :: fsvf = 89
-        INTEGER             :: i
-        
+        INTEGER(iwp)        :: fsvf = 89
+        INTEGER(iwp)        :: i
+
         DO  i = 0, io_blocks-1
             IF ( i == io_group )  THEN
@@ -4072 +4266 @@
 
                 WRITE ( fsvf )  usm_version
-                WRITE ( fsvf )  nsvfcsfl, nsvfl
+                WRITE ( fsvf )  nsvfl, ncsfl
                 WRITE ( fsvf )  svf
                 WRITE ( fsvf )  svfsurf
+                IF ( plant_canopy )  THEN
+                    WRITE ( fsvf )  csf
+                    WRITE ( fsvf )  csfsurf
+                ENDIF
                 WRITE ( fsvf )  TRIM(svf_code)
-                
+
                 CLOSE (fsvf)
 #if defined( __parallel )
-                CALL mpi_barrier( comm2d, ierr )
+                CALL MPI_BARRIER( comm2d, ierr )
 #endif
             ENDIF
@@ -4086 +4284 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Parin for &usm_par for urban surface model
-!------------------------------------------------------------------------------!
-    SUBROUTINE usm_parin
-
-       IMPLICIT NONE
-
-       CHARACTER (LEN=80) ::  line  !< string containing current line of file PARIN
-
-       NAMELIST /urban_surface_par/                                            &  
-                           land_category,                                      &
-                           mrt_factors,                                        &
-                           nrefsteps,                                          &
-                           pedestrant_category,                                &
-                           ra_horiz_coef,                                      &
-                           read_svf_on_init,                                   &
-                           roof_category,                                      &
-                           split_diffusion_radiation,                          &
-                           urban_surface,                                      &
-                           usm_anthropogenic_heat,                             &
-                           usm_energy_balance_land,                            &
-                           usm_energy_balance_wall,                            &
-                           usm_material_model,                                 &
-                           usm_lad_rma,                                        &
-                           wall_category,                                      &
-                           write_svf_on_init
-
-       line = ' '
-       
-!
-!--    Try to find urban surface model package
-       REWIND ( 11 )
-       line = ' '
-       DO   WHILE ( INDEX( line, '&urban_surface_par' ) == 0 )
-          READ ( 11, '(A)', END=10 )  line
-       ENDDO
-       BACKSPACE ( 11 )
-
-!
-!--    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
-    
-    END SUBROUTINE usm_parin    
-
-    
 !------------------------------------------------------------------------------!
 !
@@ -4146 +4290 @@
 !> Block of auxiliary subroutines:
 !> 1. quicksort and corresponding comparison
-!> 2. usm_merge_and_grow_csf for implementation of "dynamical growing" array
-!>    for svf and csf
+!> 2. usm_merge_and_grow_csf for implementation of "dynamical growing"
+!>    array for csf
 !------------------------------------------------------------------------------!    
     PURE FUNCTION svf_lt(svf1,svf2) result (res)
@@ -4168 +4312 @@
         IMPLICIT NONE
         TYPE(t_svf), DIMENSION(:), INTENT(INOUT)  :: svfl
-        INTEGER, INTENT(IN)                       :: first, last
+        INTEGER(iwp), INTENT(IN)                  :: first, last
         TYPE(t_svf)                               :: x, t
-        INTEGER                                   :: i, j
+        INTEGER(iwp)                              :: i, j
 
         IF ( first>=last ) RETURN
@@ -4217 +4361 @@
         IMPLICIT NONE
         TYPE(t_csf), DIMENSION(:), INTENT(INOUT)  :: csfl
-        INTEGER, INTENT(IN)                       :: first, last
+        INTEGER(iwp), INTENT(IN)                  :: first, last
         TYPE(t_csf)                               :: x, t
-        INTEGER                                   :: i, j
+        INTEGER(iwp)                              :: i, j
 
         IF ( first>=last ) RETURN
@@ -4243 +4387 @@
     
     SUBROUTINE usm_merge_and_grow_csf(newsize)
-        INTEGER(iwp), INTENT(in)            :: newsize  !< new array size after grow, must be >= ncsfl
-                                                        !< or -1 to shrink to minimum
-        INTEGER(iwp)                        :: iread, iwrite
-        TYPE(t_csf), DIMENSION(:), POINTER  :: acsfnew
+        INTEGER(iwp), INTENT(in)                :: newsize  !< new array size after grow, must be >= ncsfl
+                                                            !< or -1 to shrink to minimum
+        INTEGER(iwp)                            :: iread, iwrite
+        TYPE(t_csf), DIMENSION(:), POINTER      :: acsfnew
 
         IF ( newsize == -1 )  THEN
@@ -4325 +4469 @@
         INTEGER(iwp), DIMENSION(:,:), INTENT(INOUT)  :: kpcsflt
         REAL(wp), DIMENSION(:,:), INTENT(INOUT)      :: pcsflt
-        INTEGER, INTENT(IN)                          :: first, last
+        INTEGER(iwp), INTENT(IN)                     :: first, last
         REAL(wp), DIMENSION(ndcsf)                   :: t2
         INTEGER(iwp), DIMENSION(kdcsf)               :: x, t1
-        INTEGER                                      :: i, j
+        INTEGER(iwp)                                 :: i, j
 
         IF ( first>=last ) RETURN