Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3377)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3378)
@@ -28,4 +28,10 @@
 ! -----------------
 ! $Id$
+! merge from radiation branch (r3362) into trunk
+! (moh.hefny):
+! - removed read/write_svf_on_init and read_dist_max_svf (not used anymore)
+! - bugfix nzut > nzpt in calculating maxboxes
+! 
+! 3372 2018-10-18 14:03:19Z raasch
 ! bugfix: kind type of 2nd argument of mpi_win_allocate changed, misplaced
 !         __parallel directive
@@ -619,5 +625,5 @@
                 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
+                                                        !< When it switched off, only the effect of buildings and trees shadow
                                                         !< 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
@@ -776,6 +782,6 @@
                                              rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
                                              rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 
-                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m) 
-                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m) 
+                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2) 
+                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m2) 
 
     REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
@@ -987,6 +993,4 @@
     INTEGER(iwp), PARAMETER                        ::  gasize = 100000  !< initial size of growing arrays
     REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
-    REAL(wp)                                       ::  dist_max_svf = -9999.0 !< maximum distance to calculate the minimum svf to be considered. It is
-                                                                        !< used to avoid very small SVFs resulting from too far surfaces with mutual visibility
     INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
     INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
@@ -1158,6 +1162,5 @@
            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
            zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
-           write_svf_on_init, read_svf_on_init,                                &
-           nrefsteps, dist_max_svf, nsvfl, svf,                                &
+           nrefsteps, nsvfl, svf,                                              &
            svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
            surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir,      &
@@ -1167,5 +1170,5 @@
            iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
            iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
-           nsurf_type, nzub, nzut, nzpt, nzu, pch, nsurf,                      &
+           nsurf_type, nzub, nzut, nzu, pch, nsurf,                            &
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
            radiation_interactions, startwall, startland, endland, endwall,     &
@@ -2910,8 +2913,6 @@
                                   radiation_scheme, skip_time_do_radiation,    &
                                   sw_radiation, unscheduled_radiation_calls,   &
-                                  read_svf_on_init, write_svf_on_init,         &
                                   max_raytracing_dist, min_irrf_value,         &
                                   nrefsteps, raytrace_mpi_rma,                 &
-                                  dist_max_svf,                                &
                                   surface_reflections, svfnorm_report_thresh,  &
                                   radiation_interactions_on,                   &
@@ -2929,5 +2930,4 @@
                                   max_raytracing_dist, min_irrf_value,         &
                                   nrefsteps, raytrace_mpi_rma,                 &
-                                  dist_max_svf,                                &
                                   surface_reflections, svfnorm_report_thresh,  &
                                   radiation_interactions_on,                   &
@@ -4629,9 +4629,4 @@
      REAL(wp), PARAMETER               :: alpha = 0._wp      !< grid rotation (TODO: add to namelist or remove)
      REAL(wp)                          :: pc_box_area, pc_abs_frac, pc_abs_eff
-!     REAL(wp)                          ::  count_surfaces    !< number of all surfaces in model domain
-!     REAL(wp)                          ::  count_surfaces_l  !< number of all surfaces in sub-domain
-!     REAL(wp)                          ::  pt_surf_urb       !< mean surface temperature of all surfaces in model domain, temporal work-around
-!     REAL(wp)                          ::  pt_surf_urb_l     !< mean surface temperature of all surfaces in sub-domain, temporal work-around
-
      REAL(wp), DIMENSION(0:nsurf_type) :: facearea
      REAL(wp)                          :: pabsswl  = 0.0_wp  !< total absorbed SW radiation energy in local processor (W)
@@ -4686,9 +4681,10 @@
 !--         precompute effective box depth with prototype Leaf Area Density
             pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1
-            CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift),      &
+            CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift),       &
                                 60, prototype_lad,                          &
                                 CSHIFT(ABS(sunorig), pc_box_dimshift),      &
                                 pc_box_area, pc_abs_frac)
-            pc_box_area = pc_box_area * ABS(sunorig(pc_box_dimshift+1) / sunorig(1))
+            pc_box_area = pc_box_area * ABS(sunorig(pc_box_dimshift+1)      &
+                          / sunorig(1))
             pc_abs_eff = LOG(1._wp - pc_abs_frac) / prototype_lad
          ENDIF
@@ -4838,5 +4834,6 @@
            j = surfl(iy, isurf)
            i = surfl(ix, isurf)
-           surfinswdir(isurf) = rad_sw_in_dir(j,i) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
+           surfinswdir(isurf) = rad_sw_in_dir(j,i) *                        &
+                costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
         ENDDO
 !
@@ -5248,53 +5245,4 @@
       t_rad_urb = ( (pemitlw - pabslw + emissivity_urb*pinlw) / &
            (emissivity_urb*sigma_sb * area_hor) )**0.25_wp
-!
-!--  It has been turned out that the effective radiative temperature is far 
-!--  too high during nighttime, resulting in unphysical radiative forcing
-!--  with wrong signs. For the moment, as a work-around, compute the mean 
-!--  surface temperature from all surface elements, resulting in more 
-!--  physically meaningful radiative forcings.           
-!      pt_surf_urb_l    = 0.0_wp
-!      count_surfaces_l = 0.0_wp
-!      DO  m = 1, surf_lsm_h%ns
-!         k                = surf_lsm_h%k(m)
-!         pt_surf_urb_l    = pt_surf_urb_l + surf_lsm_h%pt_surface(m)            &
-!                            * exner(k)
-!         count_surfaces_l = count_surfaces_l + 1.0_wp
-!      ENDDO
-!      DO  m = 1, surf_usm_h%ns
-!         k                = surf_usm_h%k(m)
-!         pt_surf_urb_l    = pt_surf_urb_l + surf_usm_h%pt_surface(m)            &
-!                            * exner(k)
-!         count_surfaces_l = count_surfaces_l + 1.0_wp
-!      ENDDO
-!      DO  l = 0, 3
-!         DO  m = 1, surf_lsm_v(l)%ns
-!            k                = surf_lsm_v(l)%k(m)
-!            pt_surf_urb_l    = pt_surf_urb_l + surf_lsm_v(l)%pt_surface(m)      &
-!                            * exner(k)
-!            count_surfaces_l = count_surfaces_l + 1.0_wp
-!         ENDDO
-!         DO  m = 1, surf_usm_v(l)%ns
-!            k                = surf_usm_v(l)%k(m)
-!            pt_surf_urb_l    = pt_surf_urb_l + surf_usm_v(l)%pt_surface(m) * exner(k)
-!            count_surfaces_l = count_surfaces_l + 1.0_wp
-!         ENDDO
-!      ENDDO
-!      
-!      pt_surf_urb    = 0.0_wp
-!      count_surfaces = 0.0_wp
-!      
-! #if defined( __parallel )
-!      CALL MPI_ALLREDUCE( count_surfaces_l, count_surfaces, 1, MPI_REAL,        &
-!                          MPI_SUM, comm2d, ierr)
-!      CALL MPI_ALLREDUCE( pt_surf_urb_l,  pt_surf_urb,      1, MPI_REAL,        &
-!                          MPI_SUM, comm2d, ierr)
-! #else
-!      count_surfaces_l = count_surfaces
-!      pt_surf_urb_l    = pt_surf_urb
-! #endif     
-! 
-!      t_rad_urb = pt_surf_urb / count_surfaces
-     
 
     CONTAINS
@@ -5469,8 +5417,4 @@
 #endif
 
-
-       !INTEGER(iwp), DIMENSION(1:4,inorth_b:iwest_b)  ::  ijdb                               !< start and end of the local domain border coordinates (set in code)
-       !LOGICAL, DIMENSION(inorth_b:iwest_b)           ::  isborder                           !< is PE on the border of the domain in four corresponding directions
-
 !
 !--    Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be
@@ -5566,6 +5510,14 @@
 !--    check max_raytracing_dist relative to urban surface layer height
        mrl = 2.0_wp * nzu * dz(1)
+!--    set max_raytracing_dist to double the urban surface layer height, if not set
        IF ( max_raytracing_dist == -999.0_wp ) THEN
           max_raytracing_dist = mrl
+       ENDIF
+!--    check if max_raytracing_dist set too low (here we only warn the user. Other
+!      option is to correct the value again to double the urban surface layer height)
+       IF ( max_raytracing_dist  <  mrl ) THEN
+          WRITE(message_string, '(a,f6.1)') 'Max_raytracing_dist is set less than ', &
+               'double the urban surface layer height, i.e. ', mrl
+          CALL message('radiation_interaction_init', 'PA0521', 0, 0, -1, 6, 0)
        ENDIF
 !        IF ( max_raytracing_dist <= mrl ) THEN
@@ -6303,6 +6255,6 @@
                ENDDO
 
-               !
-               !--Advance itarget indices
+!
+!--            Advance itarget indices
                itarg0 = itarg1 + 1
                itarg1 = itarg1 + nzn
@@ -7467,5 +7419,5 @@
 !--      
          maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nzub,          &
-                                  nzpt - CEILING(origin(1)-.5_wp))) * nrays
+                                  nzut - CEILING(origin(1)-.5_wp))) * nrays
          IF ( ncsfl + maxboxes > ncsfla )  THEN
 !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
@@ -8284,4 +8236,7 @@
 !> It is called outside from usm_urban_surface_mod whenever the radiation fluxes
 !> are needed.
+!>
+!> This routine is not used so far. However, it may serve as an interface for radiation
+!> fluxes of urban and land surfaces
 !>
 !> TODO:
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3377)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3378)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! merge from radiation branch (r3362) into trunk
+! (moh.hefny):
+! Bugfix in the if statement to call radiation_interaction
+! 
+! 3347 2018-10-15 14:21:08Z suehring
 ! - offline nesting separated from large-scale forcing module
 ! - changes for synthetic turbulence generator
@@ -1149,5 +1154,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 3377)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3378)
@@ -28,4 +28,11 @@
 ! -----------------
 ! $Id$
+! merge from radiation branch (r3362) into trunk
+! (moh.hefny):
+! - check the requested output variables if they are correct
+! - added unscheduled_radiation_calls switch to control force_radiation_call
+! - minor formate changes
+! 
+! 3371 2018-10-18 13:40:12Z knoop
 ! Set flag indicating that albedo at urban surfaces is already initialized
 ! 
@@ -399,5 +406,5 @@
                idcsf, ndcsf, kdcsf, pct,                                       &
                startland, endland, startwall, endwall, skyvf, skyvft, nzub,    &
-               nzut, nzpt, npcbl, pcbl
+               nzut, npcbl, pcbl, unscheduled_radiation_calls
 
     USE statistics,                                                            &
@@ -2399,10 +2406,61 @@
         IMPLICIT NONE
  
-        CHARACTER (len=*),INTENT(IN)    ::  variable !:
-        CHARACTER (len=*),INTENT(OUT)   ::  unit     !:
-        
-        CHARACTER (len=varnamelength)   :: var
-
+        CHARACTER(LEN=*),INTENT(IN)    ::  variable !:
+        CHARACTER(LEN=*),INTENT(OUT)   ::  unit     !:
+
+        INTEGER(iwp)                                  :: i,j          !< index
+        CHARACTER(LEN=varnamelength)                  :: var
+        INTEGER(iwp), PARAMETER                       :: nl1 = 32     !< number of directional usm variables
+        CHARACTER(LEN=varnamelength), DIMENSION(nl1)  :: varlist1 = & !< list of directional usm variables
+                  (/'usm_rad_net','usm_rad_insw','usm_rad_inlw','usm_rad_inswdir',     &
+                    'usm_rad_inswdif','usm_rad_inswref','usm_rad_inlwdif','usm_wshf',  &
+                    'usm_rad_inlwref','usm_rad_outsw','usm_rad_outlw','usm_rad_hf',    &
+                    'usm_rad_ressw','usm_rad_reslw','usm_wghf','usm_wghf_window',      &
+                    'usm_wghf_green','usm_iwghf','usm_iwghf_window','usm_surfz',       &
+                    'usm_surfwintrans','usm_surfcat','usm_surfalb','usm_surfemis',     &
+                    'usm_t_surf','usm_t_surf_window','usm_t_surf_green','usm_t_green', &
+                    'usm_t_surf_10cm','usm_t_wall','usm_t_window','usm_t_green'/)
+
+        INTEGER(iwp), PARAMETER                       :: nl2 = 9      !< number of other variables
+        CHARACTER(LEN=varnamelength), DIMENSION(nl2)  :: varlist2 = & !< list of other usm variables
+                  (/'usm_skyvf','usm_skyvft','usm_svf','usm_dif','usm_rad_pc_inlw',    &
+                  'usm_rad_pc_insw','usm_rad_pc_inswdir','usm_rad_pc_inswdif',         &
+                  'usm_rad_pc_inswref'/)
+
+        INTEGER(iwp), PARAMETER                       :: nd = 5     !< number of directions
+        CHARACTER(LEN=6), DIMENSION(nd), PARAMETER  :: dirname = &  !< direction names
+                  (/'_roof ','_south','_north','_west ','_east '/)
+        LOGICAL                                       :: lfound     !< flag if the variable is found
+
+
+        lfound = .FALSE.
         var = TRIM(variable)
+
+!--     check if variable exists
+!       directional variables
+        DO i = 1, nl1
+           DO j = 1, nd
+              IF ( TRIM(var) == TRIM(varlist1(i))//TRIM(dirname(j)) ) THEN
+                 lfound = .TRUE.
+                 EXIT
+              ENDIF
+              IF ( lfound ) EXIT
+           ENDDO
+        ENDDO
+        IF ( lfound ) GOTO 10
+!       other variables
+        DO i = 1, nl2
+           IF ( TRIM(var) == TRIM(varlist2(i)) ) THEN
+              lfound = .TRUE.
+              EXIT
+           ENDIF
+        ENDDO
+        IF ( .NOT.  lfound ) THEN
+           unit = 'illegal'
+           RETURN
+        ENDIF
+10      CONTINUE
+
+!--     check if variable exists
         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.    &
@@ -2593,5 +2651,5 @@
 
           CASE ( 'usm_surfz' )
-!--           array of lw radiation falling to local surface after i-th reflection
+!--           array of surface height (z)
               IF ( idsint == iup_u )  THEN
                  DO  m = 1, surf_usm_h%ns
@@ -5332,10 +5390,4 @@
               ENDIF
            ENDDO
-!!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!
-!           t_window_v_p(l)%t = t_wall_v_p(l)%t
-!           surf_usm_v(l)%tt_window_m  = surf_usm_v(l)%tt_wall_m
-!           t_green_v_p(l)%t = t_wall_v_p(l)%t
-!           surf_usm_v(l)%tt_green_m  = surf_usm_v(l)%tt_wall_m
-!!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!
         ENDDO
 
@@ -6914,7 +6966,5 @@
         INTEGER(iwp)                                          :: wtc
         REAL(wp), DIMENSION(n_surface_params)                 :: wtp
-        
         LOGICAL                                               :: ascii_file = .FALSE.
-    
         INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg)   :: usm_par
         REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg)       :: usm_val
@@ -7742,7 +7792,8 @@
 !--        in case of fast changes in the skin temperature, it is required to
 !--        update the radiative fluxes in order to keep the solution stable
-           IF ( ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp ) .OR. &
+           IF ( ( ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp ) .OR.           &
                 ( ABS( t_surf_green_h_p(m) - t_surf_green_h(m) ) > 1.0_wp ) .OR. &
-                ( ABS( t_surf_window_h_p(m) - t_surf_window_h(m) ) > 1.0_wp ) ) THEN
+                ( ABS( t_surf_window_h_p(m) - t_surf_window_h(m) ) > 1.0_wp ) )  &
+                 .AND.  unscheduled_radiation_calls  ) THEN
               force_radiation_call_l = .TRUE.
            ENDIF
@@ -7897,5 +7948,5 @@
 
 !--           add RK3 term
-              t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) *       &
+              t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) *                   &
                                 surf_usm_v(l)%tt_surface_m(m)
               t_surf_window_v_p(l)%t(m) = t_surf_window_v_p(l)%t(m) + dt_3d * tsc(3) *     &
@@ -7907,6 +7958,6 @@
 !--           store also vpt_surface, which is, due to the lack of moisture on roofs simply
 !--           assumed to be the surface temperature.     
-              surf_usm_v(l)%pt_surface(m) =  ( surf_usm_v(l)%frac(ind_veg_wall,m) * t_surf_v_p(l)%t(m)  &
-                                      + surf_usm_v(l)%frac(ind_wat_win,m) * t_surf_window_v_p(l)%t(m)  &
+              surf_usm_v(l)%pt_surface(m) =  ( surf_usm_v(l)%frac(ind_veg_wall,m) * t_surf_v_p(l)%t(m)   &
+                                      + surf_usm_v(l)%frac(ind_wat_win,m) * t_surf_window_v_p(l)%t(m)    &
                                       + surf_usm_v(l)%frac(ind_pav_green,m) * t_surf_green_v_p(l)%t(m) ) &
                                       / exner(k)
@@ -7920,5 +7971,5 @@
               stend_window = ( t_surf_window_v_p(l)%t(m) - t_surf_window_v(l)%t(m) - dt_3d * tsc(3) *&
                         surf_usm_v(l)%tt_surface_window_m(m) ) / ( dt_3d  * tsc(2) )
-              stend_green = ( t_surf_green_v_p(l)%t(m) - t_surf_green_v(l)%t(m) - dt_3d * tsc(3) *&
+              stend_green = ( t_surf_green_v_p(l)%t(m) - t_surf_green_v(l)%t(m) - dt_3d * tsc(3) *   &
                         surf_usm_v(l)%tt_surface_green_m(m) ) / ( dt_3d  * tsc(2) )
 
@@ -7929,11 +7980,11 @@
                     surf_usm_v(l)%tt_surface_window_m(m) = stend_window
                     surf_usm_v(l)%tt_surface_green_m(m) = stend_green
-                 ELSEIF ( intermediate_timestep_count <                        &
+                 ELSEIF ( intermediate_timestep_count <                                 &
                           intermediate_timestep_count_max )  THEN
-                    surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend +       &
+                    surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend +                &
                                      5.3125_wp * surf_usm_v(l)%tt_surface_m(m)
-                    surf_usm_v(l)%tt_surface_green_m(m) = -9.5625_wp * stend_green +       &
+                    surf_usm_v(l)%tt_surface_green_m(m) = -9.5625_wp * stend_green +    &
                                      5.3125_wp * surf_usm_v(l)%tt_surface_green_m(m)
-                    surf_usm_v(l)%tt_surface_window_m(m) = -9.5625_wp * stend_window +       &
+                    surf_usm_v(l)%tt_surface_window_m(m) = -9.5625_wp * stend_window +  &
                                      5.3125_wp * surf_usm_v(l)%tt_surface_window_m(m)
                  ENDIF
@@ -7943,7 +7994,8 @@
 !--           update the radiative fluxes in order to keep the solution stable
 
-              IF ( ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp ) .OR. &
+              IF ( ( ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp ) .OR.            &
                    ( ABS( t_surf_green_v_p(l)%t(m) - t_surf_green_v(l)%t(m) ) > 1.0_wp ) .OR.  &
-                   ( ABS( t_surf_window_v_p(l)%t(m) - t_surf_window_v(l)%t(m) ) > 1.0_wp ) ) THEN
+                   ( ABS( t_surf_window_v_p(l)%t(m) - t_surf_window_v(l)%t(m) ) > 1.0_wp ) )   &
+                     .AND.  unscheduled_radiation_calls ) THEN
                  force_radiation_call_l = .TRUE.
               ENDIF
