Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4209)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4210)
@@ -28,4 +28,10 @@
 ! -----------------
 ! $Id$
+! - Revise steering of splitting diffuse and direct radiation
+! - Bugfixes in checks
+! - Optimize mapping of radiation components onto 2D arrays, avoid unnecessary
+!   operations
+! 
+! 4208 2019-09-02 09:01:07Z suehring
 ! Bugfix in accessing albedo_pars in the clear-sky branch (merge from branch) 
 ! 
@@ -2694,5 +2700,5 @@
 !--       Currently, 2D external radiation input is not possible in
 !--       combination with topography where average radiation is used.
-          IF ( ( rad_sw_in_f%lod == 2  .OR.  rad_sw_in_f%lod == 2  .OR.      &
+          IF ( ( rad_lw_in_f%lod == 2  .OR.  rad_sw_in_f%lod == 2  .OR.      &
                  rad_sw_in_dif_f%lod == 2  )  .AND. average_radiation )  THEN
              message_string = 'External radiation with lod = 2 is currently '//&
@@ -2704,17 +2710,17 @@
 !--       of lods divided by the number of available radiation arrays must be
 !--       1 (if all are lod = 1) or 2 (if all are lod = 2). 
-          IF ( REAL( MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
+          IF ( REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) +       &
                      MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
                      MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
                      KIND = wp ) /                                              &
-                   ( MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
+                   ( MERGE( 1.0_wp, 0.0_wp, rad_lw_in_f%from_file ) +           &
                      MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
                      MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) )       &
                      /= 1.0_wp  .AND.                                           &
-               REAL( MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
+               REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) +       &
                      MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
                      MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
                      KIND = wp ) /                                              &
-                   ( MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
+                   ( MERGE( 1.0_wp, 0.0_wp, rad_lw_in_f%from_file ) +           &
                      MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
                      MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) )       &
@@ -2808,4 +2814,6 @@
        INTEGER(iwp) ::  t   !< index of current timestep
        INTEGER(iwp) ::  tm  !< index of previous timestep 
+       
+       LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces
        
        REAL(wp) ::  fac_dt     !< interpolation factor  
@@ -2837,8 +2845,10 @@
 !--    Call clear-sky calculation for each surface orientation. 
 !--    First, horizontal surfaces
+       horizontal = .TRUE.
        surf => surf_lsm_h
        CALL radiation_external_surf
        surf => surf_usm_h
        CALL radiation_external_surf
+       horizontal = .FALSE.
 !
 !--    Vertical surfaces
@@ -3018,5 +3028,5 @@
 !--                   longwave radiation
                       IF ( ALLOCATED( rad_lw_in_diff ) )                       &
-                         rad_lw_in_diff = surf%rad_lw_in(m)
+                         rad_lw_in_diff(j,i) = surf%rad_lw_in(m)
                    ENDIF
 
@@ -3026,14 +3036,17 @@
 !
 !--          Store radiation also on 2D arrays, which are still used for 
-!--          direct-diffuse splitting.
-             DO  m = 1, surf%ns
-                i = surf%i(m)
-                j = surf%j(m)
-                
-                rad_sw_in(0,:,:)  = surf%rad_sw_in(m)
-                rad_lw_in(0,:,:)  = surf%rad_lw_in(m)
-                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
-                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
-             ENDDO
+!--          direct-diffuse splitting. Note, this is only required
+!--          for horizontal surfaces, which covers all x,y position.
+             IF ( horizontal )  THEN
+                DO  m = 1, surf%ns
+                   i = surf%i(m)
+                   j = surf%j(m)
+                   
+                   rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
+                   rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
+                   rad_sw_out(0,j,i) = surf%rad_sw_out(m)
+                   rad_lw_out(0,j,i) = surf%rad_lw_out(m)
+                ENDDO
+             ENDIF
  
           END SUBROUTINE radiation_external_surf
@@ -3052,4 +3065,7 @@
 
        INTEGER(iwp) ::  l         !< running index for surface orientation
+       
+       LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces
+       
        REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top 
        REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain 
@@ -3107,8 +3123,10 @@
 !--    Call clear-sky calculation for each surface orientation. 
 !--    First, horizontal surfaces
+       horizontal = .TRUE.
        surf => surf_lsm_h
        CALL radiation_clearsky_surf
        surf => surf_usm_h
        CALL radiation_clearsky_surf
+       horizontal = .FALSE.
 !
 !--    Vertical surfaces
@@ -3215,13 +3233,16 @@
 
 !
-!--          Fill out values in radiation arrays
-             DO  m = 1, surf%ns
-                i = surf%i(m)
-                j = surf%j(m)
-                rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
-                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
-                rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
-                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
-             ENDDO
+!--          Fill out values in radiation arrays. Note, this is only required
+!--          for horizontal surfaces, which covers all x,y position.
+             IF ( horizontal )  THEN
+                DO  m = 1, surf%ns
+                   i = surf%i(m)
+                   j = surf%j(m)
+                   rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
+                   rad_sw_out(0,j,i) = surf%rad_sw_out(m)
+                   rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
+                   rad_lw_out(0,j,i) = surf%rad_lw_out(m)
+                ENDDO
+             ENDIF
  
           END SUBROUTINE radiation_clearsky_surf
@@ -3241,4 +3262,6 @@
 
        INTEGER(iwp) ::  l         !< running index for surface orientation
+       
+       LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces
 
        REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top 
@@ -3285,8 +3308,10 @@
 !
 !--    First, horizontal surfaces
+       horizontal = .TRUE.
        surf => surf_lsm_h
        CALL radiation_constant_surf
        surf => surf_usm_h
        CALL radiation_constant_surf
+       horizontal = .FALSE.
 !
 !--    Vertical surfaces
@@ -3397,13 +3422,16 @@
 
 !
-!--          Fill out values in radiation arrays
-             DO  m = 1, surf%ns
-                i = surf%i(m)
-                j = surf%j(m)
-                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
-                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
-                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
-                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
-             ENDDO
+!--          Fill out values in radiation arrays. Note, this is only required
+!--          for horizontal surfaces, which covers all x,y position.
+             IF ( horizontal )  THEN
+                DO  m = 1, surf%ns
+                   i = surf%i(m)
+                   j = surf%j(m)
+                   rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
+                   rad_sw_out(0,j,i) = surf%rad_sw_out(m)
+                   rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
+                   rad_lw_out(0,j,i) = surf%rad_lw_out(m)
+                ENDDO
+             ENDIF
 
           END SUBROUTINE radiation_constant_surf
@@ -5421,7 +5449,15 @@
 !--  Split downwelling shortwave radiation into a diffuse and a direct part. 
 !--  Note, if radiation scheme is RRTMG or diffuse radiation is externally
-!--  prescribed, this is not required. 
-     IF (  radiation_scheme /= 'rrtmg'  .AND.                                  &
-           .NOT. rad_sw_in_dif_f%from_file )  CALL calc_diffusion_radiation
+!--  prescribed, this is not required. Please note, in case of external
+!--  radiation, the clear-sky model is applied during spinup, so that 
+!--  radiation need to be split also in this case. 
+     IF ( radiation_scheme == 'constant'   .OR.                                &
+          radiation_scheme == 'clear-sky'  .OR.                                &
+          ( radiation_scheme == 'external'  .AND.                              &
+            .NOT. rad_sw_in_dif_f%from_file  )  .OR.                           &
+          ( radiation_scheme == 'external'  .AND.                              &
+            time_since_reference_point < 0.0_wp ) )  THEN
+        CALL calc_diffusion_radiation
+     ENDIF
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -6110,5 +6146,5 @@
       area_hor = pinlw / rad_lw_in_diff(nyn,nxl)
       t_rad_urb = ( ( pemitlw - pabslw + emissivity_urb * pinlw ) / &
-           (emissivity_urb * sigma_sb * area_hor) )**0.25_wp
+           (emissivity_urb * sigma_sb * area_hor) )**0.25_wp     
 
      IF ( debug_output_timestep )  CALL debug_message( 'radiation_interaction', 'end' )
