Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 2975)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 2977)
@@ -25,4 +25,11 @@
 ! -----------------
 ! $Id$
+! Implement changes from branch radiation (r2948-2971) with minor modifications
+! (moh.hefny):
+! - set radiation_interactions according to the existence of urban/land vertical
+!   surfaces and trees to activiate RTM
+! - set average_radiation to TRUE if RTM is activiated
+! 
+! 2938 2018-03-27 15:52:42Z suehring
 ! - Revise Inifor initialization for geostrophic wind components
 ! - Initialize synthetic turbulence generator in case of Inifor initialization  
@@ -506,9 +513,10 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  radiation_init, radiation, radiation_control, radiation_scheme, &
+        ONLY:  average_radiation,                                              &
+               radiation_init, radiation, radiation_control, radiation_scheme, &
                radiation_calc_svf, radiation_write_svf,                        &
                radiation_interaction, radiation_interactions,                  &
                radiation_interaction_init, radiation_read_svf,                 &
-               radiation_presimulate_solar_pos
+               radiation_presimulate_solar_pos, radiation_interactions_on
     
     USE random_function_mod 
@@ -534,5 +542,5 @@
     USE surface_mod,                                                           &
         ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
-                surf_usm_h, get_topography_top_index_ji
+                surf_usm_h, get_topography_top_index_ji, vertical_surfaces_exist
    
     USE transpose_indices
@@ -2357,5 +2365,5 @@
 !-- If required, set chemical emissions
 !-- (todo(FK): This should later on be CALLed time-dependently in init_3d_model)
-    IF ( air_chemistry ) THEN
+    IF ( air_chemistry )  THEN
        CALL chem_emissions
     ENDIF
@@ -2364,4 +2372,23 @@
 !-- Initialize radiation processes
     IF ( radiation )  THEN
+!
+!--    Activate radiation_interactions according to the existence of vertical surfaces and/or trees.
+!--    The namelist parameter radiation_interactions_on can override this behavior.
+!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in 
+!--    init_surface_arrays.)
+       IF ( radiation_interactions_on )  THEN
+          IF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
+             radiation_interactions    = .TRUE.
+             average_radiation         = .TRUE.
+          ELSE
+             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions 
+                                                   !< calculations necessary in case of flat surface
+          ENDIF
+       ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
+          message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
+                           'vertical surfaces and/or trees exist. The model will run ' // &
+                           'without RTM (no shadows, no radiation reflections)'
+          CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
+       ENDIF
 !
 !--    If required, initialize radiation interactions between surfaces 
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 2975)
+++ palm/trunk/SOURCE/palm.f90	(revision 2977)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Deduct spinup_time from RUN_CONTROL output of main 3d run 
+! (use time_since_reference_point instead of simulated_time)
+! 
+! 2951 2018-04-06 09:05:08Z kanani
 ! Add log_point_s for pmci_init
 ! 
@@ -240,4 +244,5 @@
                nest_domain, neutral, nudging, passive_scalar, runnr,           &
                simulated_time, simulated_time_chr, spinup,                     &
+               time_since_reference_point,                                     &
                user_interface_current_revision,                                &
                user_interface_required_revision, version, wall_heatflux,       &
@@ -464,5 +469,5 @@
 !
 !-- Set start time in format hh:mm:ss
-    simulated_time_chr = time_to_string( simulated_time )
+    simulated_time_chr = time_to_string( time_since_reference_point )
 
 !
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2975)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2977)
@@ -19,5 +19,5 @@
 !
 ! Current revisions:
-! -----------------
+! ------------------
 ! 
 ! 
@@ -25,4 +25,11 @@
 ! -----------------
 ! $Id$
+! Implement changes from branch radiation (r2948-2971) with minor modifications,
+! plus some formatting.
+! (moh.hefny):
+! Add plant canopy type to account for changes in LAD (based on the changes
+! done by Resler & Pavel) and correct the error message to PALM Standard.
+! 
+! 2932 2018-03-26 09:39:22Z maronga
 ! renamed canopy_par to plant_canopy_parameters
 ! 
@@ -213,4 +220,7 @@
     REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
     REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
+
+    REAL(wp) ::  lad_type_coef(0:10) = 1.0_wp   !< multiplicative coeficients for particular types
+                                                !< of plant canopy (e.g. deciduous tree during winter)
 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
@@ -962,5 +972,5 @@
                                   alpha_lad, beta_lad, canopy_drag_coeff,      &
                                   canopy_mode, cthf,                           &
-                                  lad_surface,                                 & 
+                                  lad_surface, lad_type_coef,                  & 
                                   lad_vertical_gradient,                       &
                                   lad_vertical_gradient_level,                 &
@@ -971,5 +981,5 @@
        NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
                                   canopy_mode, cthf,                           &
-                                  lad_surface,                                 & 
+                                  lad_surface, lad_type_coef,                  & 
                                   lad_vertical_gradient,                       &
                                   lad_vertical_gradient_level,                 &
@@ -1034,7 +1044,7 @@
 !>
 !> num_levels
-!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
-!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
-!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
+!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
+!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
+!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
 !> ...
 !>
@@ -1060,4 +1070,5 @@
 
        INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
+       INTEGER(iwp)                        ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
        INTEGER(iwp)                        ::  i, j      !< running index
        INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
@@ -1073,29 +1084,35 @@
 !
 !--    Open and read plant canopy input data
-       OPEN(152, file='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
-                 access='SEQUENTIAL', action='READ', status='OLD',             &
-                 form='FORMATTED', err=515)
-       READ(152, *, err=516, end=517) nzp   !< read first line = number of vertical layers
+       OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
+                 ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD',             &
+                 FORM='FORMATTED', ERR=515)
+       READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
        
-       ALLOCATE(col(0:nzp-1))
+       ALLOCATE( col(0:nzp-1) )
 
        DO
-          READ(152, *, err=516, end=517) dtype, i, j, col(:)
-          IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
-
-             SELECT CASE (dtype)
-                CASE( 1 )   !< leaf area density
-!
-!--                 This is just the pure canopy layer assumed to be grounded to
-!--                 a flat domain surface. At locations where plant canopy sits 
-!--                 on top of any kind of topography, the vertical plant column 
-!--                 must be "lifted", which is done in SUBROUTINE pcm_tendency.
-                    lad_s(0:nzp-1, j, i) = col(0:nzp-1)
-                    
-                CASE DEFAULT
-                   WRITE(message_string, '(a,i2,a)')   &
-                      'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
-                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
-             END SELECT
+          READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
+          IF ( i < nxlg  .OR.  i > nxrg  .OR.  j < nysg  .OR.  j > nyng )  CYCLE
+          
+          SELECT CASE (dtype)
+             CASE( 1 )   !< leaf area density
+!
+!--             This is just the pure canopy layer assumed to be grounded to
+!--             a flat domain surface. At locations where plant canopy sits 
+!--             on top of any kind of topography, the vertical plant column 
+!--             must be "lifted", which is done in SUBROUTINE pcm_tendency.           
+                IF ( pctype < 0  .OR.  pctype > 10 )  THEN   !< incorrect plant canopy type
+                   WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //  &
+                                              'Allowed values 0 <= pctype <= 10, ' //  &
+                                              'but pctype is ', pctype
+                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
+                ENDIF
+                lad_s(0:nzp-1,j,i) = col(0:nzp-1) * lad_type_coef(pctype)
+                
+             CASE DEFAULT
+                WRITE(message_string, '(a,i2,a)')   &
+                     'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
+                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
+          END SELECT
        ENDDO
 
@@ -1107,5 +1124,5 @@
 
 517    CLOSE(152)
-       DEALLOCATE(col)
+       DEALLOCATE( col )
        
        CALL exchange_horiz( lad_s, nbgp )
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2975)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2977)
@@ -22,5 +22,5 @@
 !
 ! Current revisions:
-! -----------------
+! ------------------
 ! 
 ! 
@@ -28,4 +28,19 @@
 ! -----------------
 ! $Id$
+! Implement changes from branch radiation (r2948-2971) with minor modifications,
+! plus some formatting.
+! (moh.hefny):
+! - replaced plant_canopy by npcbl to check tree existence to avoid weird 
+!   allocation of related arrays (after domain decomposition some domains 
+!   contains no trees although plant_canopy (global parameter) is still TRUE).
+! - added a namelist parameter to force RTM settings
+! - enabled the option to switch radiation reflections off
+! - renamed surf_reflections to surface_reflections 
+! - removed average_radiation flag from the namelist (now it is implicitly set 
+!   in init_3d_model according to RTM)
+! - edited read and write sky view factors and CSF routines to account for
+!   the sub-domains which may not contain any of them
+! 
+! 2967 2018-04-13 11:22:08Z raasch
 ! bugfix: missing parallel cpp-directives added
 ! 
@@ -432,8 +447,9 @@
                 sun_direction = .FALSE.,              & !< flag parameter indicating whether solar direction shall be calculated
                 average_radiation = .FALSE.,          & !< flag to set the calculation of radiation averaging for the domain
-                radiation_interactions = .FALSE.,     & !< flag to control if radiation interactions via sky-view factors shall be considered
-                surf_reflections = .TRUE.               !< flag to switch the calculation of radiation interaction between surfaces.
+                radiation_interactions = .FALSE.,     & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
+                surface_reflections = .TRUE.,         & !< flag to switch the calculation of radiation interaction between surfaces.
                                                         !< When it switched off, only the effect of buildings and trees shadow will
                                                         !< will be considered. However fewer SVFs are expected.
+                radiation_interactions_on = .TRUE.      !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
 
 
@@ -680,5 +696,5 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pct              !< top layer of the plant canopy
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pch              !< heights of the plant canopy
-    INTEGER(iwp)                                   ::  npcbl            !< number of the plant canopy gridboxes in local processor
+    INTEGER(iwp)                                   ::  npcbl = 0        !< number of the plant canopy gridboxes in local processor
     INTEGER(wp), DIMENSION(:,:), ALLOCATABLE       ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
@@ -955,5 +971,5 @@
            surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart,     &
            surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw,      &
-           pcbl, npcbl, iup_u, inorth_u, isouth_u, ieast_u, iwest_u,           &
+           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,           &
            iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
            nsurf_type, nzub, nzut, nzu, pch, nsurf,                            &
@@ -961,5 +977,5 @@
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
            radiation_interactions, startwall, startland, endland, endwall,     &
-           skyvf, skyvft
+           skyvf, skyvft, radiation_interactions_on, average_radiation
 
 #if defined ( __rrtmg )
@@ -1309,14 +1325,4 @@
              CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
           ENDIF
-       ENDIF
-
-!
-!--    Radiation interactions
-       IF ( nrefsteps < 1  .AND.  radiation_interactions )  THEN
-          message_string = 'nrefsteps must be > 0 when using LSM/USM to' //    &
-                           'account for surface outgoing SW flux.'       //    &
-                           'You may set surf_reflections = .FALSE. to '  //    &
-                           'diable surface reflections instead.'
-          CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 )
        ENDIF
 
@@ -2634,6 +2640,6 @@
                                   nrefsteps, mrt_factors, rma_lad_raytrace,    &
                                   dist_max_svf,                                &
-                                  average_radiation,                           &
-                                  surf_reflections, svfnorm_report_thresh
+                                  surface_reflections, svfnorm_report_thresh,  &
+                                  radiation_interactions_on
    
        NAMELIST /radiation_parameters/   albedo, albedo_type, albedo_lw_dir,   &
@@ -2647,6 +2653,6 @@
                                   nrefsteps, mrt_factors, rma_lad_raytrace,    &
                                   dist_max_svf,                                &
-                                  average_radiation,                           &
-                                  surf_reflections, svfnorm_report_thresh
+                                  surface_reflections, svfnorm_report_thresh,  &
+                                  radiation_interactions_on
    
        line = ' '
@@ -2693,10 +2699,11 @@
        radiation = .TRUE.
 
+       IF ( .NOT.  radiation_interactions_on  .AND.  surface_reflections )  THEN
+          message_string = 'surface_reflections is allowed only when '      // &
+               'radiation_interactions_on is set to TRUE'
+          CALL message( 'radiation_parin', 'PA0487',1, 2, 0, 6, 0 )
+       ENDIF
+
  12    CONTINUE
-       
-
-
-!--    Set radiation_interactions flag according to urban_ and land_surface flag
-       IF ( urban_surface  .OR.  land_surface )  radiation_interactions = .TRUE.
        
     END SUBROUTINE radiation_parin
@@ -4354,5 +4361,5 @@
          sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
 
-         IF ( plant_canopy )  THEN
+         IF ( npcbl > 0 )  THEN
 !--            precompute effective box depth with prototype Leaf Area Density
             pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1
@@ -4465,14 +4472,16 @@
 #endif
 
-     DO isvf = 1, nsvfl
-         isurf = svfsurf(1, isvf)
-         k = surfl(iz, isurf)
-         j = surfl(iy, isurf)
-         i = surfl(ix, isurf)
-         isurfsrc = svfsurf(2, isvf)
-
-!--         for surface-to-surface factors we calculate thermal radiation in 1st pass
-         surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
-     ENDDO
+     IF ( surface_reflections)  THEN
+        DO  isvf = 1, nsvfl
+           isurf = svfsurf(1, isvf)
+           k     = surfl(iz, isurf)
+           j     = surfl(iy, isurf)
+           i     = surfl(ix, isurf)
+           isurfsrc = svfsurf(2, isvf)
+!
+!--        For surface-to-surface factors we calculate thermal radiation in 1st pass
+           surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
+        ENDDO
+     ENDIF
 
      !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors
@@ -4494,5 +4503,5 @@
      ENDIF
 
-     IF ( plant_canopy )  THEN
+     IF ( npcbl > 0 )  THEN
 
          pcbinswdir(:) = 0._wp
@@ -4536,4 +4545,14 @@
 !        surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw
 
+     IF ( .NOT.  surface_reflections )  THEN
+!
+!--     Set nrefsteps to 0 to disable reflections        
+        nrefsteps = 0
+        surfoutsl = albedo_surf * surfins
+        surfoutll = (1._wp - emiss_surf) * surfinl
+        surfoutsw = surfoutsw + surfoutsl
+        surfoutlw = surfoutlw + surfoutll
+     ENDIF
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !--     Next passes - reflections
@@ -4585,5 +4604,5 @@
 
 !--  push heat flux absorbed by plant canopy to respective 3D arrays
-     IF ( plant_canopy )  THEN
+     IF ( npcbl > 0 )  THEN
          pc_heating_rate(:,:,:) = 0._wp
          DO ipcgb = 1, npcbl
@@ -4733,74 +4752,71 @@
 !--  domain when using average_radiation in the respective radiation model
 
-     IF ( average_radiation )  THEN
-
-!--     precalculate face areas for all face directions using normal vector
-        DO d = 0, nsurf_type
-            facearea(d) = 1._wp
-            IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
-            IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
-            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
-        ENDDO
-!
-!--     absorbed/received SW & LW and emitted LW energy of all physical
-!--     surfaces (land and urban) in local processor
-        pinswl = 0._wp
-        pinlwl = 0._wp
-        pabsswl = 0._wp
-        pabslwl = 0._wp
-        pemitlwl = 0._wp
-        emiss_sum_surfl = 0._wp
-        area_surfl = 0._wp
-        DO  i = 1, nsurfl
-           d = surfl(id, i)
-!--        received SW & LW
-           pinswl = pinswl + surfinsw(i) * facearea(d)
-           pinlwl = pinlwl + surfinlw(i) * facearea(d)
-!--        absorbed SW & LW
-           pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
-                                                  surfinsw(i) * facearea(d)
-           pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
-!--        emitted LW
-           pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
-!--        emissivity and area sum
-           emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
-           area_surfl = area_surfl + facearea(d)
-        END DO
-!
-!--     add the absorbed SW energy by plant canopy
-        IF ( plant_canopy )  THEN
-           pabsswl = pabsswl + SUM(pcbinsw)
-           pabslwl = pabslwl + SUM(pcbinlw)
-        ENDIF
-!
-!--     gather all rad flux energy in all processors
+!--  Precalculate face areas for all face directions using normal vector
+     DO d = 0, nsurf_type
+        facearea(d) = 1._wp
+        IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
+        IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
+        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
+     ENDDO
+!
+!--  absorbed/received SW & LW and emitted LW energy of all physical
+!--  surfaces (land and urban) in local processor
+     pinswl = 0._wp
+     pinlwl = 0._wp
+     pabsswl = 0._wp
+     pabslwl = 0._wp
+     pemitlwl = 0._wp
+     emiss_sum_surfl = 0._wp
+     area_surfl = 0._wp
+     DO  i = 1, nsurfl
+        d = surfl(id, i)
+!--  received SW & LW
+        pinswl = pinswl + surfinsw(i) * facearea(d)
+        pinlwl = pinlwl + surfinlw(i) * facearea(d)
+!--   absorbed SW & LW
+        pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
+                                                surfinsw(i) * facearea(d)
+        pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
+!--   emitted LW
+        pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
+!--   emissivity and area sum
+        emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
+        area_surfl = area_surfl + facearea(d)
+     END DO
+!
+!--  add the absorbed SW energy by plant canopy
+     IF ( npcbl > 0 )  THEN
+        pabsswl = pabsswl + SUM(pcbinsw)
+        pabslwl = pabslwl + SUM(pcbinlw)
+     ENDIF
+!
+!--  gather all rad flux energy in all processors
 #if defined( __parallel )
-        CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
-        CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
-        CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
-        CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
-        CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
-        CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
-        CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
+     CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
+     CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
 #else
-        pinsw = pinswl
-        pinlw = pinlwl
-        pabssw = pabsswl
-        pabslwl = pabslw
-        pemitlwl = pemitlw
-        emiss_sum_surf = emiss_sum_surfl
-        area_surf = area_surfl
+     pinsw = pinswl
+     pinlw = pinlwl
+     pabssw = pabsswl
+     pabslwl = pabslw
+     pemitlwl = pemitlw
+     emiss_sum_surf = emiss_sum_surfl
+     area_surf = area_surfl
 #endif
 
-!--     (1) albedo
-        IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
-
-!--     (2) average emmsivity
-        IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf
-
-!--     (3) temperature
-        t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
-
-     ENDIF
+!--  (1) albedo
+     IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
+
+!--  (2) average emmsivity
+     IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf
+
+!--  (3) temperature
+     t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
+     
 
     CONTAINS
@@ -5071,5 +5087,5 @@
 
 !--    fill gridpcbl and pcbl
-       IF ( plant_canopy )  THEN
+       IF ( npcbl > 0 )  THEN
            ALLOCATE( pcbl(iz:ix, 1:npcbl) )
            ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) )
@@ -5593,54 +5609,57 @@
 
             DEALLOCATE ( zdirs, zbdry, vffrac, ztransp )
-
-            DO isurfs = 1, nsurf
-                IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
-                    surfl(iz, isurflt), surfl(id, isurflt), &
-                    surf(ix, isurfs), surf(iy, isurfs), &
-                    surf(iz, isurfs), surf(id, isurfs)) )  THEN
-                    CYCLE
-                ENDIF
+!
+!--         Following calculations only required for surface_reflections
+            IF ( surface_reflections )  THEN
+
+               DO  isurfs = 1, nsurf
+                  IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
+                     surfl(iz, isurflt), surfl(id, isurflt), &
+                     surf(ix, isurfs), surf(iy, isurfs), &
+                     surf(iz, isurfs), surf(id, isurfs)) )  THEN
+                     CYCLE
+                  ENDIF
                   
-                sd = surf(id, isurfs)
-                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 = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
-                sqdist = SUM(uv(:)**2)
-                uv = uv / SQRT(sqdist)
-
-!--             reject raytracing above max distance
-                IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
-                    ray_skip_maxdist = ray_skip_maxdist + 1
-                    CYCLE
-                ENDIF
-                
-!--             irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
-                rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
-                    * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
-                    / (pi * sqdist) & ! square of distance between centers
-                    * facearea(sd)
-
-!--             reject raytracing for potentially too small view factor values
-                IF ( rirrf < min_irrf_value ) THEN
-                    ray_skip_minval = ray_skip_minval + 1
-                    CYCLE
-                ENDIF
-
-!--             raytrace + process plant canopy sinks within
-                CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
-                        visible, transparency, win_lad)
-
-                IF ( .NOT.  visible ) CYCLE
-                ! rsvf = rirrf * transparency
-
-!--             write to the svf array
-                nsvfl = nsvfl + 1
-!--             check dimmension of asvf array and enlarge it if needed
-                IF ( nsvfla < nsvfl )  THEN
-                    k = nsvfla * 2
-                    IF ( msvf == 0 )  THEN
+                  sd = surf(id, isurfs)
+                  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 = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
+                  sqdist = SUM(uv(:)**2)
+                  uv = uv / SQRT(sqdist)
+
+!--               reject raytracing above max distance
+                  IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
+                     ray_skip_maxdist = ray_skip_maxdist + 1
+                     CYCLE
+                  ENDIF
+                 
+!--               irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
+                  rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
+                      * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
+                      / (pi * sqdist) & ! square of distance between centers
+                      * facearea(sd)
+
+!--               reject raytracing for potentially too small view factor values
+                  IF ( rirrf < min_irrf_value ) THEN
+                      ray_skip_minval = ray_skip_minval + 1
+                      CYCLE
+                  ENDIF
+
+!--               raytrace + process plant canopy sinks within
+                  CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
+                                visible, transparency, win_lad)
+
+                  IF ( .NOT.  visible ) CYCLE
+                 ! rsvf = rirrf * transparency
+
+!--               write to the svf array
+                  nsvfl = nsvfl + 1
+!--               check dimmension of asvf array and enlarge it if needed
+                  IF ( nsvfla < nsvfl )  THEN
+                     k = nsvfla * 2
+                     IF ( msvf == 0 )  THEN
                         msvf = 1
                         ALLOCATE( asvf1(k) )
@@ -5648,5 +5667,5 @@
                         asvf1(1:nsvfla) = asvf2
                         DEALLOCATE( asvf2 )
-                    ELSE
+                     ELSE
                         msvf = 0
                         ALLOCATE( asvf2(k) )
@@ -5654,17 +5673,18 @@
                         asvf2(1:nsvfla) = asvf1
                         DEALLOCATE( asvf1 )
-                    ENDIF
-
-!                     WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
-!                     CALL radiation_write_debug_log( msg )
-                    
-                    nsvfla = k
-                ENDIF
-!--             write svf values into the array
-                asvf(nsvfl)%isurflt = isurflt
-                asvf(nsvfl)%isurfs = isurfs
-                asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency
-                asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor
-            ENDDO
+                     ENDIF
+
+!                      WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
+!                      CALL radiation_write_debug_log( msg )
+                     
+                     nsvfla = k
+                  ENDIF
+!--               write svf values into the array
+                  asvf(nsvfl)%isurflt = isurflt
+                  asvf(nsvfl)%isurfs = isurfs
+                  asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency
+                  asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor
+               ENDDO
+            ENDIF
         ENDDO
 
@@ -6733,5 +6753,5 @@
  !--             read nsvfl, ncsfl
               READ ( fsvf ) nsvfl, ncsfl, nsurfl_from_file
-              IF ( nsvfl <= 0  .OR.  ncsfl < 0 )  THEN
+              IF ( nsvfl < 0  .OR.  ncsfl < 0 )  THEN
                   WRITE( message_string, * ) 'Wrong number of SVF or CSF'
                   CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 )
@@ -6747,32 +6767,37 @@
                   CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 )
               ENDIF
-              
-              IF ( .NOT. ALLOCATED( skyvf ) )    ALLOCATE( skyvf(nsurfl) )
-              IF ( .NOT. ALLOCATED( skyvft ) )   ALLOCATE( skyvft(nsurfl) )
-              IF ( .NOT. ALLOCATED( svf ) )      ALLOCATE( svf(ndsvf,nsvfl) )
-              IF ( .NOT. ALLOCATED( svfsurf ) )  ALLOCATE( svfsurf(idsvf,nsvfl) )
-              READ(fsvf) skyvf
-              READ(fsvf) skyvft
-              READ(fsvf) svf
-              READ(fsvf) svfsurf
-              IF ( plant_canopy )  THEN
-                  IF ( .NOT. ALLOCATED( csf ) )      ALLOCATE( csf(ndcsf,ncsfl) )
-                  IF ( .NOT. ALLOCATED( csfsurf ) )  ALLOCATE( csfsurf(idcsf,ncsfl) )
-                  READ(fsvf) csf
-                  READ(fsvf) csfsurf
+              IF ( nsurfl > 0 )  THEN
+                 IF ( .NOT. ALLOCATED( skyvf ) )    ALLOCATE( skyvf(nsurfl) )
+                 IF ( .NOT. ALLOCATED( skyvft ) )   ALLOCATE( skyvft(nsurfl) )
+                 READ(fsvf) skyvf
+                 READ(fsvf) skyvft
+              ENDIF
+                
+              IF ( nsvfl > 0 )  THEN
+                 IF ( .NOT. ALLOCATED( svf ) )      ALLOCATE( svf(ndsvf,nsvfl) )
+                 IF ( .NOT. ALLOCATED( svfsurf ) )  ALLOCATE( svfsurf(idsvf,nsvfl) )
+                 READ(fsvf) svf
+                 READ(fsvf) svfsurf
+              ENDIF
+
+              IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
+                 IF ( .NOT. ALLOCATED( csf ) )      ALLOCATE( csf(ndcsf,ncsfl) )
+                 IF ( .NOT. ALLOCATED( csfsurf ) )  ALLOCATE( csfsurf(idcsf,ncsfl) )
+                 READ(fsvf) csf
+                 READ(fsvf) csfsurf
               ENDIF
               READ ( fsvf ) svf_code_field
-              
+                
               IF ( TRIM(svf_code_field) /= TRIM(svf_code) )  THEN
-                  WRITE( message_string, * ) 'Wrong structure of binary svf file'
-                  CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 )
+                 WRITE( message_string, * ) 'Wrong structure of binary svf file'
+                 CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 )
               ENDIF       
 !
-!--          Close binary file                 
-             CALL close_file( fsvf )
+!--           Close binary file                 
+              CALL close_file( fsvf )
                 
            ENDIF
 #if defined( __parallel )
-            CALL MPI_BARRIER( comm2d, ierr )
+           CALL MPI_BARRIER( comm2d, ierr )
 #endif
         ENDDO
@@ -6801,9 +6826,13 @@
                 WRITE ( fsvf )  rad_version
                 WRITE ( fsvf )  nsvfl, ncsfl, nsurfl
-                WRITE ( fsvf )  skyvf
-                WRITE ( fsvf )  skyvft
-                WRITE ( fsvf )  svf
-                WRITE ( fsvf )  svfsurf
-                IF ( plant_canopy )  THEN
+                IF ( nsurfl > 0 ) THEN
+                   WRITE ( fsvf )  skyvf
+                   WRITE ( fsvf )  skyvft
+                ENDIF
+                IF ( nsvfl > 0 ) THEN
+                   WRITE ( fsvf )  svf
+                   WRITE ( fsvf )  svfsurf
+                ENDIF
+                IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
                     WRITE ( fsvf )  csf
                     WRITE ( fsvf )  csfsurf
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 2975)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 2977)
@@ -26,4 +26,11 @@
 ! -----------------
 ! $Id$
+! Implement changes from branch radiation (r2948-2971) with minor modifications,
+! plus some formatting.
+! (moh.hefny)
+! Added flag to check the existence of vertical urban/land surfaces, required
+! to activate RTM
+! 
+! 2970 2018-04-13 15:09:23Z suehring
 ! Remove un-necessary initialization of surface elements in old large-scale 
 ! forcing mode
@@ -450,4 +457,7 @@
     INTEGER(iwp) ::  ns_v_on_file(0:3)                       !< total number of vertical surfaces with the same facing, required for writing restart data 
 
+    LOGICAL ::  vertical_surfaces_exist = .FALSE.   !< flag indicating that there are vertical urban/land surfaces 
+                                                    !< in the domain (required to activiate RTM)
+
 
     SAVE
@@ -488,5 +498,5 @@
     PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file,       &
            ns_v_on_file, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,       &
-           surf_usm_h, surf_usm_v, surf_type
+           surf_usm_h, surf_usm_v, surf_type, vertical_surfaces_exist
 !
 !-- Public subroutines and functions
@@ -615,4 +625,8 @@
     SUBROUTINE init_surface_arrays
 
+
+       USE pegrid
+
+
        IMPLICIT NONE
 
@@ -628,4 +642,7 @@
        INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces 
        INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces 
+
+       INTEGER(iwp)              ::  num_surf_v_l !< number of vertically-aligned local urban/land surfaces
+       INTEGER(iwp)              ::  num_surf_v   !< number of vertically-aligned total urban/land surfaces
 
        LOGICAL ::  building                       !< flag indicating building grid point
@@ -920,5 +937,18 @@
                                                nys, nyn, nxl, nxr )
        ENDDO
-
+!
+!--    Set the flag for the existence of vertical urban/land surfaces
+       num_surf_v_l = 0
+       DO  l = 0, 3
+          num_surf_v_l = num_surf_v_l + surf_usm_v(l)%ns + surf_lsm_v(l)%ns
+       ENDDO
+
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( num_surf_v_l, num_surf_v, 1, MPI_INTEGER, MPI_SUM, comm2d, ierr)
+#else
+       num_surf_v = num_surf_v_l
+#endif
+       IF ( num_surf_v > 0 ) vertical_surfaces_exist = .TRUE.
+        
 
     END SUBROUTINE init_surface_arrays
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 2975)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 2977)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Implement changes from branch radiation (r2948-2971) with minor modifications.
+! (moh.hefny)
+! Fixed bug in if statement
+! 
+! 2941 2018-04-03 11:54:58Z kanani
 ! Deduct spinup_time from RUN_CONTROL output of main 3d run 
 ! (use time_since_reference_point instead of simulated_time)
@@ -1034,5 +1039,5 @@
                 CALL cpu_log( log_point(50), 'radiation', 'stop' )
 
-                IF ( urban_surface  .OR.  land_surface  .AND.                  &
+                IF ( ( urban_surface  .OR.  land_surface )  .AND.               &
                      radiation_interactions )  THEN
                    CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2975)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2977)
@@ -28,4 +28,9 @@
 ! -----------------
 ! $Id$
+! Implement changes from branch radiation (r2948-2971) with minor modifications.
+! (moh.hefny):
+! Extended exn for all model domain height to avoid the need to get nzut.
+! 
+! 2963 2018-04-12 14:47:44Z suehring
 ! Introduce index for vegetation/wall, pavement/green-wall and water/window 
 ! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
@@ -3303,5 +3308,4 @@
             DO  i = nxl, nxr
                 DO  k = nzb_do, nzt_do
-!                   print*, j,nys,nyn,i,nxl,nxr,k,nzb_do,nzt_do,local_pf(i,j,k),temp_pf(k,j,i)
                     local_pf(i,j,k) = temp_pf(k,j,i)
                 ENDDO
