Index: palm/trunk/LIB/rrtmg/rrtmg_sw_rad.nomcica.f90
===================================================================
--- palm/trunk/LIB/rrtmg/rrtmg_sw_rad.nomcica.f90	(revision 3271)
+++ palm/trunk/LIB/rrtmg/rrtmg_sw_rad.nomcica.f90	(revision 3272)
@@ -78,13 +78,15 @@
       subroutine rrtmg_sw &
             (ncol    ,nlay    ,icld    ,iaer    , &
-             play    ,plev    ,tlay    ,tlev    ,tsfc    , &
-             h2ovmr  ,o3vmr   ,co2vmr  ,ch4vmr  ,n2ovmr  ,o2vmr, &
-             asdir   ,asdif   ,aldir   ,aldif   , &
-             coszen  ,adjes   ,dyofyr  ,scon    , &
-             inflgsw ,iceflgsw,liqflgsw,cldfr   , &
-             taucld  ,ssacld  ,asmcld  ,fsfcld  , &
-             cicewp  ,cliqwp  ,reice   ,reliq   , &
-             tauaer  ,ssaaer  ,asmaer  ,ecaer   , &
-             swuflx  ,swdflx  ,swhr    ,swuflxc ,swdflxc ,swhrc)
+             play    ,plev    ,tlay    ,tlev    , &
+             tsfc    ,h2ovmr  ,o3vmr   ,co2vmr  , &
+             ch4vmr  ,n2ovmr  ,o2vmr   ,asdir   , &
+             asdif   ,aldir   ,aldif   ,coszen  , &
+             adjes   ,dyofyr  ,scon    ,inflgsw , &
+             iceflgsw,liqflgsw,cldfr   ,taucld  , &
+             ssacld  ,asmcld  ,fsfcld  ,cicewp  , &
+             cliqwp  ,reice   ,reliq   ,tauaer  , &
+             ssaaer  ,asmaer  ,ecaer   ,swuflx  , &
+             swdflx  ,swhr    ,swuflxc ,swdflxc , &
+             swhrc   ,dirdflux,difdflux)
 
 ! ------- Description -------
@@ -169,4 +171,6 @@
 !   delta scaling based on setting of idelm flag. 
 !     Dec 2008: M. J. Iacono, AER, Inc.
+!-- Added the output direct and diffuse fluxes to the argument of rrtmg_sw
+!     Aug 2018: M. Salim, Humboldt Uni, Berlin, Germany
 
 ! --------- Modules ---------
@@ -417,6 +421,6 @@
       real(kind=rb) :: swnflx(nlay+2)         ! Total sky shortwave net flux (W/m2)
       real(kind=rb) :: swnflxc(nlay+2)        ! Clear sky shortwave net flux (W/m2)
-      real(kind=rb) :: dirdflux(nlay+2)       ! Direct downward shortwave surface flux
-      real(kind=rb) :: difdflux(nlay+2)       ! Diffuse downward shortwave surface flux
+      real(kind=rb), intent(out) :: dirdflux(nlay+2)       ! Direct downward shortwave surface flux
+      real(kind=rb), intent(out) :: difdflux(nlay+2)       ! Diffuse downward shortwave surface flux
       real(kind=rb) :: uvdflx(nlay+2)         ! Total sky downward shortwave flux, UV/vis  
       real(kind=rb) :: nidflx(nlay+2)         ! Total sky downward shortwave flux, near-IR 
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3271)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3272)
@@ -28,4 +28,17 @@
 ! -----------------
 ! $Id$
+! - split direct and diffusion shortwave radiation using RRTMG rather than using
+!   calc_diffusion_radiation, in case of RRTMG
+! - removed the namelist variable split_diffusion_radiation. Now splitting depends
+!   on the choise of radiation radiation scheme
+! - removed calculating the rdiation flux for surfaces at the radiation scheme
+!   in case of using RTM since it will be calculated anyway in the radiation 
+!   interaction routine.
+! - set SW radiation flux for surfaces to zero at night in case of no RTM is used
+! - precalculate the unit vector yxdir of ray direction to avoid the temporarly 
+!   array allocation during the subroutine call
+! - fixed a bug in calculating the max number of boxes ray can cross in the domain
+! 
+! 3264 2018-09-20 13:54:11Z moh.hefny
 ! Bugfix in raytrace_2d calls
 ! 
@@ -722,6 +735,7 @@
                                              rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
                                              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_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²) 
 
     REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
@@ -822,6 +836,4 @@
 
 !-- configuration parameters (they can be setup in PALM config)
-    LOGICAL                                        ::  split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation
-                                                                                          !< (.F. in case the radiation model already does it)    
     LOGICAL                                        ::  rma_lad_raytrace = .FALSE.         !< use MPI RMA to access LAD for raytracing (instead of global array)
     LOGICAL                                        ::  mrt_factors = .FALSE.              !< whether to generate MRT factor files during init
@@ -1080,5 +1092,4 @@
            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
            zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
-           split_diffusion_radiation,                                          &
            nrefsteps, mrt_factors, dist_max_svf, nsvfl, svf,                   &
            svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
@@ -1153,5 +1164,5 @@
        SELECT CASE ( TRIM( var ) )
 
-         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
+         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' )
              IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
                 message_string = '"output of "' // TRIM( var ) // '" requi' // &
@@ -1189,4 +1200,5 @@
              IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
              IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
              IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
              IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
@@ -2769,5 +2781,4 @@
                                   radiation_scheme, skip_time_do_radiation,    &
                                   sw_radiation, unscheduled_radiation_calls,   &
-                                  split_diffusion_radiation,                   &
                                   max_raytracing_dist, min_irrf_value,         &
                                   nrefsteps, mrt_factors, rma_lad_raytrace,    &
@@ -2782,5 +2793,4 @@
                                   radiation_scheme, skip_time_do_radiation,    &
                                   sw_radiation, unscheduled_radiation_calls,   &
-                                  split_diffusion_radiation,                   &
                                   max_raytracing_dist, min_irrf_value,         &
                                   nrefsteps, mrt_factors, rma_lad_raytrace,    &
@@ -3054,5 +3064,5 @@
                 rad_lw_out(k,:,:) = rrtm_lwuflx(0,k)
              ENDDO
-
+             rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
 !
 !--          Save heating rates (convert from K/d to K/h). 
@@ -3070,55 +3080,30 @@
              ENDDO
 
-!
-!--          Save surface radiative fluxes and change in LW heating rate 
-!--          onto respective surface elements
-!--          Horizontal surfaces
-             IF ( surf_lsm_h%ns > 0 )  THEN
-                surf_lsm_h%rad_lw_in           = rrtm_lwdflx(0,nzb)
-                surf_lsm_h%rad_lw_out          = rrtm_lwuflx(0,nzb)
-                surf_lsm_h%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
-             ENDIF             
-             IF ( surf_usm_h%ns > 0 )  THEN
-                surf_usm_h%rad_lw_in           = rrtm_lwdflx(0,nzb)
-                surf_usm_h%rad_lw_out          = rrtm_lwuflx(0,nzb)
-                surf_usm_h%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
-             ENDIF
-!
-!--          Vertical surfaces. 
-             DO  l = 0, 3
-                IF ( surf_lsm_v(l)%ns > 0 )  THEN
-                   surf_lsm_v(l)%rad_lw_in           = rrtm_lwdflx(0,nzb)
-                   surf_lsm_v(l)%rad_lw_out          = rrtm_lwuflx(0,nzb)
-                   surf_lsm_v(l)%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
-                ENDIF
-                IF ( surf_usm_v(l)%ns > 0 )  THEN
-                   surf_usm_v(l)%rad_lw_in           = rrtm_lwdflx(0,nzb)
-                   surf_usm_v(l)%rad_lw_out          = rrtm_lwuflx(0,nzb)
-                   surf_usm_v(l)%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)
-                ENDIF
-             ENDDO
-
           ENDIF
 
           IF ( sw_radiation .AND. sun_up )  THEN
-             CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer        ,&
-             rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
-             rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
-             rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir       ,&
-             rrtm_asdif      , rrtm_aldir   , rrtm_aldif , zenith,           &
-             0.0_wp          , day_of_year  , solar_constant,   rrtm_inflgsw,&
-             rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
-             rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
-             rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
-             rrtm_sw_ssaaer  , rrtm_sw_asmaer  , rrtm_sw_ecaer ,             &
-             rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
-             rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
+             CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld     , rrtm_iaer      ,&
+             rrtm_play      , rrtm_plev     , rrtm_tlay     , rrtm_tlev      ,&
+             rrtm_tsfc      , rrtm_h2ovmr   , rrtm_o3vmr    , rrtm_co2vmr    ,&
+             rrtm_ch4vmr    , rrtm_n2ovmr   , rrtm_o2vmr    , rrtm_asdir     ,&
+             rrtm_asdif     , rrtm_aldir    , rrtm_aldif    , zenith         ,&
+             0.0_wp         , day_of_year   , solar_constant, rrtm_inflgsw   ,&
+             rrtm_iceflgsw  , rrtm_liqflgsw , rrtm_cldfr    , rrtm_sw_taucld ,&
+             rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp    ,&
+             rrtm_cliqwp    , rrtm_reice    , rrtm_reliq    , rrtm_sw_tauaer ,&
+             rrtm_sw_ssaaer , rrtm_sw_asmaer, rrtm_sw_ecaer , rrtm_swuflx    ,&
+             rrtm_swdflx    , rrtm_swhr     , rrtm_swuflxc  , rrtm_swdflxc   ,&
+             rrtm_swhrc     , rrtm_dirdflux , rrtm_difdflux )
   
 !
-!--          Save fluxes
+!--          Save fluxes:
+!--          - whole domain
              DO k = nzb, nzt+1
                 rad_sw_in(k,:,:)  = rrtm_swdflx(0,k)
                 rad_sw_out(k,:,:) = rrtm_swuflx(0,k)
              ENDDO
+!--          - direct and diffuse SW at urban-surface-layer (required by RTM)
+             rad_sw_in_dir(:,:) = rrtm_dirdflux(0,nzb)
+             rad_sw_in_diff(:,:) = rrtm_difdflux(0,nzb)
 
 !
@@ -3128,29 +3113,4 @@
                 rad_sw_cs_hr(k,:,:)  = rrtm_swhrc(0,k) * d_hours_day
              ENDDO
-
-!
-!--          Save surface radiative fluxes onto respective surface elements
-!--          Horizontal surfaces
-             IF ( surf_lsm_h%ns > 0 )  THEN
-                   surf_lsm_h%rad_sw_in     = rrtm_swdflx(0,nzb)
-                   surf_lsm_h%rad_sw_out    = rrtm_swuflx(0,nzb)
-             ENDIF
-             IF ( surf_usm_h%ns > 0 )  THEN
-                   surf_usm_h%rad_sw_in     = rrtm_swdflx(0,nzb)
-                   surf_usm_h%rad_sw_out    = rrtm_swuflx(0,nzb)
-             ENDIF
-!
-!--          Vertical surfaces. Fluxes are obtain at respective vertical 
-!--          level of the surface element
-             DO  l = 0, 3
-                IF ( surf_lsm_v(l)%ns > 0 )  THEN
-                      surf_lsm_v(l)%rad_sw_in  = rrtm_swdflx(0,nzb)
-                      surf_lsm_v(l)%rad_sw_out = rrtm_swuflx(0,nzb)
-                ENDIF             
-                IF ( surf_usm_v(l)%ns > 0 )  THEN
-                      surf_usm_v(l)%rad_sw_in  = rrtm_swdflx(0,nzb)
-                      surf_usm_v(l)%rad_sw_out = rrtm_swuflx(0,nzb)
-                ENDIF       
-             ENDDO
 !
 !--       Solar radiation is zero during night
@@ -3158,8 +3118,10 @@
              rad_sw_in  = 0.0_wp
              rad_sw_out = 0.0_wp
+             rad_sw_in_dir(:,:) = 0.0_wp
+             rad_sw_in_diff(:,:) = 0.0_wp
           ENDIF
 !
 !--    RRTMG is called for each (j,i) grid point separately, starting at the 
-!--    highest topography level
+!--    highest topography level. Here no RTM is used since average_radiation is false
        ELSE
 !
@@ -3551,5 +3513,7 @@
                                   rrtm_swuflxc(:,k_topo:nzt_rad+1),            & 
                                   rrtm_swdflxc(:,k_topo:nzt_rad+1),            &
-                                  rrtm_swhrc(:,k_topo+1:nzt_rad+1) )
+                                  rrtm_swhrc(:,k_topo+1:nzt_rad+1),            &
+                                  rrtm_dirdflux(:,k_topo:nzt_rad+1),           &
+                                  rrtm_difdflux(:,k_topo:nzt_rad+1) )
 
                    DEALLOCATE( rrtm_sw_taucld_dum )
@@ -3609,4 +3573,35 @@
                    rad_sw_in  = 0.0_wp
                    rad_sw_out = 0.0_wp
+!--             !!!!!!!! ATTENSION !!!!!!!!!!!!!!!
+!--             Surface radiative fluxes should be also set to zero here                  
+!--                Save surface radiative fluxes onto respective surface elements
+!--                Horizontal surfaces
+                   DO  m = surf_lsm_h%start_index(j,i),                        &
+                           surf_lsm_h%end_index(j,i)
+                      surf_lsm_h%rad_sw_in(m)     = 0.0_wp
+                      surf_lsm_h%rad_sw_out(m)    = 0.0_wp
+                   ENDDO             
+                   DO  m = surf_usm_h%start_index(j,i),                        &
+                           surf_usm_h%end_index(j,i)
+                      surf_usm_h%rad_sw_in(m)     = 0.0_wp
+                      surf_usm_h%rad_sw_out(m)    = 0.0_wp
+                   ENDDO 
+!
+!--                Vertical surfaces. Fluxes are obtain at respective vertical 
+!--                level of the surface element
+                   DO  l = 0, 3
+                      DO  m = surf_lsm_v(l)%start_index(j,i),                  &
+                              surf_lsm_v(l)%end_index(j,i)
+                         k                           = surf_lsm_v(l)%k(m)
+                         surf_lsm_v(l)%rad_sw_in(m)  = 0.0_wp
+                         surf_lsm_v(l)%rad_sw_out(m) = 0.0_wp
+                      ENDDO             
+                      DO  m = surf_usm_v(l)%start_index(j,i),                  &
+                              surf_usm_v(l)%end_index(j,i)
+                         k                           = surf_usm_v(l)%k(m)
+                         surf_usm_v(l)%rad_sw_in(m)  = 0.0_wp
+                         surf_usm_v(l)%rad_sw_out(m) = 0.0_wp
+                      ENDDO 
+                   ENDDO
                 ENDIF
 
@@ -3617,34 +3612,36 @@
 !
 !--    Finally, calculate surface net radiation for surface elements.
-!--    First, for horizontal surfaces    
-       DO  m = 1, surf_lsm_h%ns
-          surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m)                      &
-                                - surf_lsm_h%rad_sw_out(m)                     &
-                                + surf_lsm_h%rad_lw_in(m)                      &
-                                - surf_lsm_h%rad_lw_out(m)
-       ENDDO
-       DO  m = 1, surf_usm_h%ns
-          surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m)                      &
-                                - surf_usm_h%rad_sw_out(m)                     &
-                                + surf_usm_h%rad_lw_in(m)                      &
-                                - surf_usm_h%rad_lw_out(m)
-       ENDDO
-!
-!--    Vertical surfaces. 
-!--    Todo: weight with azimuth and zenith angle according to their orientation!
-       DO  l = 0, 3     
-          DO  m = 1, surf_lsm_v(l)%ns
-             surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m)             &
-                                      - surf_lsm_v(l)%rad_sw_out(m)            &
-                                      + surf_lsm_v(l)%rad_lw_in(m)             &
-                                      - surf_lsm_v(l)%rad_lw_out(m)
+       IF (  .NOT.  radiation_interactions  ) THEN
+!--       First, for horizontal surfaces    
+          DO  m = 1, surf_lsm_h%ns
+             surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m)                   &
+                                   - surf_lsm_h%rad_sw_out(m)                  &
+                                   + surf_lsm_h%rad_lw_in(m)                   &
+                                   - surf_lsm_h%rad_lw_out(m)
           ENDDO
-          DO  m = 1, surf_usm_v(l)%ns
-             surf_usm_v(l)%rad_net(m) = surf_usm_v(l)%rad_sw_in(m)             &
-                                      - surf_usm_v(l)%rad_sw_out(m)            &
-                                      + surf_usm_v(l)%rad_lw_in(m)             &
-                                      - surf_usm_v(l)%rad_lw_out(m)
+          DO  m = 1, surf_usm_h%ns
+             surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m)                   &
+                                   - surf_usm_h%rad_sw_out(m)                  &
+                                   + surf_usm_h%rad_lw_in(m)                   &
+                                   - surf_usm_h%rad_lw_out(m)
           ENDDO
-       ENDDO
+!
+!--       Vertical surfaces. 
+!--       Todo: weight with azimuth and zenith angle according to their orientation!
+          DO  l = 0, 3     
+             DO  m = 1, surf_lsm_v(l)%ns
+                surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m)          &
+                                         - surf_lsm_v(l)%rad_sw_out(m)         &
+                                         + surf_lsm_v(l)%rad_lw_in(m)          &
+                                         - surf_lsm_v(l)%rad_lw_out(m)
+             ENDDO
+             DO  m = 1, surf_usm_v(l)%ns
+                surf_usm_v(l)%rad_net(m) = surf_usm_v(l)%rad_sw_in(m)          &
+                                         - surf_usm_v(l)%rad_sw_out(m)         &
+                                         + surf_usm_v(l)%rad_lw_in(m)          &
+                                         - surf_usm_v(l)%rad_lw_out(m)
+             ENDDO
+          ENDDO
+       ENDIF
 
 
@@ -3905,4 +3902,6 @@
           DEALLOCATE ( rrtm_swhr  )
           DEALLOCATE ( rrtm_swhrc )
+          DEALLOCATE ( rrtm_dirdflux )
+          DEALLOCATE ( rrtm_difdflux )
 
        ENDIF
@@ -4070,4 +4069,6 @@
        ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
        ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_dirdflux(0:0,nzb:nzt_rad+1) )
+       ALLOCATE ( rrtm_difdflux(0:0,nzb:nzt_rad+1) )
 
        rrtm_swdflx  = 0.0_wp
@@ -4077,4 +4078,6 @@
        rrtm_swdflxc = 0.0_wp
        rrtm_swhrc   = 0.0_wp
+       rrtm_dirdflux = 0.0_wp
+       rrtm_difdflux = 0.0_wp
 
        ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
@@ -4482,5 +4485,6 @@
 
      IMPLICIT NONE
-     INTEGER(iwp)                      :: i, j, k, kk, d, refstep, m, mm, l, ll
+
+     INTEGER(iwp)                      :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
      INTEGER(iwp)                      :: isurf, isurfsrc, isvf, icsf, ipcgb
      INTEGER(iwp)                      :: isd                !< solar direction number
@@ -4519,5 +4523,4 @@
      REAL(wp)                          :: area_hor           !< total horizontal area of domain in all processor
 
-
 #if ! defined( __nopointer )
      IF ( plant_canopy )  THEN
@@ -4545,13 +4548,13 @@
 
      IF ( zenith(0) > 0 )  THEN
-!--         now we will "squash" the sunorig vector by grid box size in
-!--         each dimension, so that this new direction vector will allow us
-!--         to traverse the ray path within grid coordinates directly
+!--      now we will "squash" the sunorig vector by grid box size in
+!--      each dimension, so that this new direction vector will allow us
+!--      to traverse the ray path within grid coordinates directly
          sunorig_grid = (/ sunorig(1)/dz(1), sunorig(2)/dy, sunorig(3)/dx /)
-!--         sunorig_grid = sunorig_grid / norm2(sunorig_grid)
+!--      sunorig_grid = sunorig_grid / norm2(sunorig_grid)
          sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
 
          IF ( npcbl > 0 )  THEN
-!--            precompute effective box depth with prototype Leaf Area Density
+!--         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),      &
@@ -4564,14 +4567,7 @@
      ENDIF
 
-!--     split diffusion and direct part of the solar downward radiation
-!--     comming from radiation model and store it in 2D arrays
-!--     rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff
-     IF ( split_diffusion_radiation )  THEN
-         CALL calc_diffusion_radiation
-     ELSE
-         rad_sw_in_diff = 0.0_wp
-         rad_sw_in_dir(:,:)  = rad_sw_in(0,:,:)
-         rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
-     ENDIF
+!--  if radiation scheme is not RRTMG, split diffusion and direct part of the solar downward radiation
+!--  comming from radiation model and store it in 2D arrays
+     IF (  radiation_scheme /= 'rrtmg'  )  CALL calc_diffusion_radiation
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -6696,12 +6692,12 @@
          !--
          maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1)
-!          IF ( ncsfl + maxboxes > ncsfla )  THEN
-! !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
-! !--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
-! !--                                                / log(grow_factor)), kind=wp))
-! !--         or use this code to simply always keep some extra space after growing
-!             k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
-!             CALL merge_and_grow_csf(k)
-!          ENDIF
+         IF ( ncsfl + maxboxes > ncsfla )  THEN
+!--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
+!--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
+!--                                                / log(grow_factor)), kind=wp))
+!--         or use this code to simply always keep some extra space after growing
+            k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
+            CALL merge_and_grow_csf(k)
+         ENDIF
 
          !--Calculate transparencies and store new CSFs
@@ -6748,7 +6744,4 @@
                      IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
                         curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
-                        IF (ncsfl+1 > SIZE(acsf) )                             &
-                           CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
-
                         ncsfl = ncsfl + 1
                         acsf(ncsfl)%ip = ip
@@ -6810,6 +6803,4 @@
 
                      IF ( create_csf )  THEN
-                        IF (ncsfl+1 > SIZE(acsf) )                             &
-                           CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
                         ncsfl = ncsfl + 1
                         acsf(ncsfl)%ip = ip
