Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4189)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4190)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! type real_1d changed to real_1d_3d
+! 
+! 4186 2019-08-23 16:06:14Z suehring
 ! Minor formatting adjustments
 ! 
@@ -396,12 +399,15 @@
     END TYPE int_2d_32bit
 !
-!-- Define data type to read 1D real variables
-    TYPE real_1d
+!-- Define data type to read 1D or 3D real variables. 
+    TYPE real_1d_3d
        LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
 
+       INTEGER(iwp) ::  lod = -1        !< level-of-detail
+       
        REAL(wp) ::  fill = -9999.9_wp                  !< fill value
        
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  var     !< respective variable
-    END TYPE real_1d   
+       REAL(wp), DIMENSION(:),     ALLOCATABLE ::  var1d     !< respective 1D variable
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var3d     !< respective 3D variable
+    END TYPE real_1d_3d   
 !
 !-- Define data type to read 2D real variables
@@ -409,4 +415,6 @@
        LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
 
+       INTEGER(iwp) ::  lod             !< level-of-detail
+       
        REAL(wp) ::  fill = -9999.9_wp                !< fill value
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var !< respective variable
@@ -747,10 +755,12 @@
 !
 !-- Public data structures
-    PUBLIC real_1d,                                                            &
+    PUBLIC real_1d_3d,                                                         &
            real_2d
 !
 !-- Public variables
     PUBLIC albedo_pars_f, albedo_type_f, basal_area_density_f, buildings_f,    &
-           building_id_f, building_pars_f, building_type_f, char_fill,         &
+           building_id_f, building_pars_f, building_type_f,                    &
+           char_fill,                                                          &
+           char_lod,                                                           &
            chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
            coord_ref_sys,                                                      &
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4189)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4190)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Implement external radiation forcing also for level-of-detail = 2 
+! (horizontally 2D radiation)
+! 
+! 4188 2019-08-26 14:15:47Z suehring
 ! Minor adjustment in error message
 ! 
@@ -235,4 +239,5 @@
                water_type_f,                                                   &
                char_fill,                                                      &
+               char_lod,                                                       &
                check_existence,                                                &
                close_input_file,                                               &
@@ -247,5 +252,5 @@
                pids_id,                                                        &
                open_read_file,                                                 &
-               real_1d,                                                        &
+               real_1d_3d,                                                     &
                vars_pids
 
@@ -780,10 +785,10 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
 !
-!-- External radiation
-    TYPE( real_1d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model 
-    TYPE( real_1d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model 
-    TYPE( real_1d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
-    TYPE( real_1d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model 
-
+!-- External radiation. Depending on the given level of detail either a 1D or 
+!-- a 3D array will be allocated.
+    TYPE( real_1d_3d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model 
+    TYPE( real_1d_3d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model 
+    TYPE( real_1d_3d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
+    TYPE( real_1d_3d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model 
 
     INTERFACE radiation_check_data_output
@@ -1458,5 +1463,5 @@
        INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces 
        INTEGER(iwp) ::  m         !< running index for surface elements
-       INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps 
+       INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps
 #if defined( __rrtmg )
        INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
@@ -2501,10 +2506,11 @@
           IF ( check_existence( vars_pids, 'time_rad' ) )  THEN
              CALL netcdf_data_input_get_dimension_length( pids_id,             &
-                                                          ntime, 'time_rad' )
+                                                          ntime,               &
+                                                          'time_rad' )
           
-             ALLOCATE( time_rad_f%var(0:ntime-1) )
+             ALLOCATE( time_rad_f%var1d(0:ntime-1) )
 !                                                                                 
 !--          Read variable                    
-             CALL get_variable( pids_id, 'time_rad', time_rad_f%var )
+             CALL get_variable( pids_id, 'time_rad', time_rad_f%var1d )
                                 
              time_rad_f%from_file = .TRUE.
@@ -2513,42 +2519,95 @@
 !--       Input shortwave downwelling.
           IF ( check_existence( vars_pids, 'rad_sw_in' ) )  THEN
-             ALLOCATE( rad_sw_in_f%var(0:ntime-1) )
 !         
-!--          Read _FillValue attribute
+!--          Get _FillValue attribute
              CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill,         &
                                  .FALSE., 'rad_sw_in' )
-!                                                                                 
-!--          Read variable                    
-             CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var )
-                                
-             rad_sw_in_f%from_file = .TRUE.
+!         
+!--          Get level-of-detail
+             CALL get_attribute( pids_id, char_lod, rad_sw_in_f%lod,           &
+                                 .FALSE., 'rad_sw_in' )
+!
+!--          Level-of-detail 1 - radiation depends only on time_rad
+             IF ( rad_sw_in_f%lod == 1 )  THEN
+                ALLOCATE( rad_sw_in_f%var1d(0:ntime-1) )
+                CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var1d )
+                rad_sw_in_f%from_file = .TRUE.
+!
+!--          Level-of-detail 2 - radiation depends on time_rad, y, x
+             ELSEIF ( rad_sw_in_f%lod == 2 )  THEN  
+                ALLOCATE( rad_sw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
+                
+                CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var3d,    &
+                                   nxl, nxr, nys, nyn, 0, ntime-1 )
+                                   
+                rad_sw_in_f%from_file = .TRUE.
+             ELSE 
+                message_string = '"rad_sw_in" has no valid lod attribute'
+                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
+             ENDIF
           ENDIF
 !         
 !--       Input longwave downwelling.
           IF ( check_existence( vars_pids, 'rad_lw_in' ) )  THEN
-             ALLOCATE( rad_lw_in_f%var(0:ntime-1) )
 !         
-!--          Read _FillValue attribute
+!--          Get _FillValue attribute
              CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill,         &
                                  .FALSE., 'rad_lw_in' )
-!                                                                                 
-!--          Read variable                    
-             CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var )
-                                
-             rad_lw_in_f%from_file = .TRUE.
+!         
+!--          Get level-of-detail
+             CALL get_attribute( pids_id, char_lod, rad_lw_in_f%lod,           &
+                                 .FALSE., 'rad_lw_in' )
+!
+!--          Level-of-detail 1 - radiation depends only on time_rad
+             IF ( rad_lw_in_f%lod == 1 )  THEN
+                ALLOCATE( rad_lw_in_f%var1d(0:ntime-1) )
+                CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var1d )
+                rad_lw_in_f%from_file = .TRUE.
+!
+!--          Level-of-detail 2 - radiation depends on time_rad, y, x
+             ELSEIF ( rad_lw_in_f%lod == 2 )  THEN  
+                ALLOCATE( rad_lw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
+                
+                CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var3d,    &
+                                   nxl, nxr, nys, nyn, 0, ntime-1 )
+                                   
+                rad_lw_in_f%from_file = .TRUE.
+             ELSE 
+                message_string = '"rad_lw_in" has no valid lod attribute'
+                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
+             ENDIF
           ENDIF
 !         
 !--       Input shortwave downwelling, diffuse part.
           IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) )  THEN
-             ALLOCATE( rad_sw_in_dif_f%var(0:ntime-1) )
 !         
 !--          Read _FillValue attribute
              CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill,     &
                                  .FALSE., 'rad_sw_in_dif' )
-!                                                                                 
-!--          Read variable                    
-             CALL get_variable( pids_id, 'rad_sw_in_dif', rad_sw_in_dif_f%var )
-                                
-             rad_sw_in_dif_f%from_file = .TRUE.
+!         
+!--          Get level-of-detail
+             CALL get_attribute( pids_id, char_lod, rad_sw_in_dif_f%lod,       &
+                                 .FALSE., 'rad_sw_in_dif' )
+!
+!--          Level-of-detail 1 - radiation depends only on time_rad
+             IF ( rad_sw_in_dif_f%lod == 1 )  THEN
+                ALLOCATE( rad_sw_in_dif_f%var1d(0:ntime-1) )
+                CALL get_variable( pids_id, 'rad_sw_in_dif',                   &
+                                   rad_sw_in_dif_f%var1d )
+                rad_sw_in_dif_f%from_file = .TRUE.
+!
+!--          Level-of-detail 2 - radiation depends on time_rad, y, x
+             ELSEIF ( rad_sw_in_dif_f%lod == 2 )  THEN  
+                ALLOCATE( rad_sw_in_dif_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
+                
+                CALL get_variable( pids_id, 'rad_sw_in_dif',                   &
+                                   rad_sw_in_dif_f%var3d,                      &
+                                   nxl, nxr, nys, nyn, 0, ntime-1 )
+                                   
+                rad_sw_in_dif_f%from_file = .TRUE.
+             ELSE 
+                message_string = '"rad_sw_in_dif" has no valid lod attribute'
+                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
+             ENDIF
           ENDIF
 !         
@@ -2573,26 +2632,5 @@
           ENDIF
           
-          IF ( ANY( rad_sw_in_f%var ==  rad_sw_in_f%fill ) )  THEN
-             message_string = 'External radiation forcing: rad_sw_in must ' // &
-                              'not contain any fill values'
-             CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
-          ENDIF
-          
-          IF ( ANY( rad_lw_in_f%var ==  rad_lw_in_f%fill ) )  THEN
-             message_string = 'External radiation forcing: rad_lw_in must ' // &
-                              'not contain any fill values'
-             CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
-          ENDIF
-          
-          IF ( rad_sw_in_dif_f%from_file )  THEN
-             IF ( ANY( rad_sw_in_dif_f%var ==  rad_sw_in_dif_f%fill ) )  THEN
-                message_string = 'External radiation forcing: ' //             &
-                                 'rad_lw_in_dif must not contain any fill ' // &
-                                 'values'
-                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
-             ENDIF
-          ENDIF
-          
-          IF ( time_rad_f%var(0) /= 0.0_wp )  THEN
+          IF ( time_rad_f%var1d(0) /= 0.0_wp )  THEN
              message_string = 'External radiation forcing: first point in ' // &
                               'time is /= 0.0.'
@@ -2600,10 +2638,92 @@
           ENDIF
           
-          IF ( end_time - spinup_time > time_rad_f%var(ntime-1) )  THEN
+          IF ( end_time - spinup_time > time_rad_f%var1d(ntime-1) )  THEN
              message_string = 'External radiation forcing does not cover ' //  &
                               'the entire simulation time.'
              CALL message( 'radiation_init', 'PA0314', 1, 2, 0, 6, 0 )
           ENDIF
+!
+!--       Check for fill values in radiation
+          IF ( ALLOCATED( rad_sw_in_f%var1d ) )  THEN
+             IF ( ANY( rad_sw_in_f%var1d == rad_sw_in_f%fill ) )  THEN
+                message_string = 'External radiation array "rad_sw_in" ' //    &
+                                 'must not contain any fill values.'
+                CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
           
+          IF ( ALLOCATED( rad_lw_in_f%var1d ) )  THEN
+             IF ( ANY( rad_lw_in_f%var1d == rad_lw_in_f%fill ) )  THEN
+                message_string = 'External radiation array "rad_lw_in" ' //    &
+                                 'must not contain any fill values.'
+                CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+          
+          IF ( ALLOCATED( rad_sw_in_dif_f%var1d ) )  THEN
+             IF ( ANY( rad_sw_in_dif_f%var1d == rad_sw_in_dif_f%fill ) )  THEN
+                message_string = 'External radiation array "rad_sw_in_dif" ' //&
+                                 'must not contain any fill values.'
+                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+          
+          IF ( ALLOCATED( rad_sw_in_f%var3d ) )  THEN
+             IF ( ANY( rad_sw_in_f%var3d == rad_sw_in_f%fill ) )  THEN
+                message_string = 'External radiation array "rad_sw_in" ' //    &
+                                 'must not contain any fill values.'
+                CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+          
+          IF ( ALLOCATED( rad_lw_in_f%var3d ) )  THEN
+             IF ( ANY( rad_lw_in_f%var3d == rad_lw_in_f%fill ) )  THEN
+                message_string = 'External radiation array "rad_lw_in" ' //    &
+                                 'must not contain any fill values.'
+                CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+          
+          IF ( ALLOCATED( rad_sw_in_dif_f%var3d ) )  THEN
+             IF ( ANY( rad_sw_in_dif_f%var3d == rad_sw_in_dif_f%fill ) )  THEN
+                message_string = 'External radiation array "rad_sw_in_dif" ' //&
+                                 'must not contain any fill values.'
+                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+!
+!--       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.      &
+                 rad_sw_in_dif_f%lod == 2  )  .AND. average_radiation )  THEN
+             message_string = 'External radiation with lod = 2 is currently '//&
+                              'not possible with average_radiation = .T..'
+                CALL message( 'radiation_init', 'PA0670', 1, 2, 0, 6, 0 )
+          ENDIF
+!
+!--       All radiation input should have the same level of detail. The sum 
+!--       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 ) +       &
+                     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_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 ) +       &
+                     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_sw_in_f%from_file ) +           &
+                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) )       &
+                     /= 2.0_wp )  THEN
+             message_string = 'External radiation input should have the same '//&
+                              'lod.'
+             CALL message( 'radiation_init', 'PA0673', 1, 2, 0, 6, 0 )
+          ENDIF
+
        ENDIF
 !
@@ -2690,8 +2810,5 @@
        INTEGER(iwp) ::  tm  !< index of previous timestep 
        
-       REAL(wp) ::  fac_dt     !< interpolation factor
-       REAL(wp) ::  lw_in      !< downwelling longwave radiation, interpolated value     
-       REAL(wp) ::  sw_in      !< downwelling shortwave radiation, interpolated value       
-       REAL(wp) ::  sw_in_dif  !< downwelling diffuse shortwave radiation, interpolated value   
+       REAL(wp) ::  fac_dt     !< interpolation factor  
 
        TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine   
@@ -2708,5 +2825,5 @@
        ELSE
           t = 0
-          DO WHILE ( time_rad_f%var(t) <= time_since_reference_point )
+          DO WHILE ( time_rad_f%var1d(t) <= time_since_reference_point )
              t = t + 1
           ENDDO
@@ -2714,19 +2831,8 @@
           tm = MAX( t-1, 0 )
           
-          fac_dt = ( time_since_reference_point - time_rad_f%var(tm) + dt_3d ) &
-                 / ( time_rad_f%var(t)  - time_rad_f%var(tm) )
+          fac_dt = ( time_since_reference_point - time_rad_f%var1d(tm) + dt_3d ) &
+                 / ( time_rad_f%var1d(t)  - time_rad_f%var1d(tm) )
           fac_dt = MIN( 1.0_wp, fac_dt )
        ENDIF
-             
-       sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var(tm)                       &
-                          + fac_dt * rad_sw_in_f%var(t)
-                                         
-       lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var(tm)                       &
-                          + fac_dt * rad_lw_in_f%var(t)
-!
-!--    Limit shortwave incoming radiation to positive values, in order 
-!--    to overcome possible observation errors.
-       sw_in = MAX( 0.0_wp, sw_in )
-       sw_in = MERGE( sw_in, 0.0_wp, sun_up )
 !
 !--    Call clear-sky calculation for each surface orientation. 
@@ -2744,21 +2850,4 @@
           CALL radiation_external_surf
        ENDDO
-!
-!--    In case of average radiation check if external diffuse shortwave 
-!--    radiation is available and, if required, interpolate it onto the
-!--    respective arrays.
-       IF ( average_radiation )  THEN
-          IF ( rad_sw_in_dif_f%from_file )  THEN 
-             sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var(tm)          &
-                                   + fac_dt * rad_sw_in_dif_f%var(t)
-                                         
-             IF ( ALLOCATED( rad_sw_in_diff ) )  rad_sw_in_diff = sw_in_dif
-             IF ( ALLOCATED( rad_sw_in_dir  ) )  rad_sw_in_dir  = sw_in - sw_in_dif
-!
-!--          Diffuse longwave radiation equals the total downwelling longwaver
-!--          radiation
-             IF ( ALLOCATED( rad_lw_in_diff ) )  rad_lw_in_diff = lw_in
-          ENDIF
-       ENDIF
        
        CONTAINS
@@ -2770,29 +2859,96 @@
              IMPLICIT NONE
 
-             INTEGER(iwp) ::  i         !< grid index along x-dimension   
-             INTEGER(iwp) ::  j         !< grid index along y-dimension  
-             INTEGER(iwp) ::  k         !< grid index along z-dimension    
-             INTEGER(iwp) ::  m         !< running index for surface elements       
+             INTEGER(iwp) ::  i    !< grid index along x-dimension   
+             INTEGER(iwp) ::  j    !< grid index along y-dimension  
+             INTEGER(iwp) ::  k    !< grid index along z-dimension    
+             INTEGER(iwp) ::  m    !< running index for surface elements
+             
+             REAL(wp) ::  lw_in     !< downwelling longwave radiation, interpolated value     
+             REAL(wp) ::  sw_in     !< downwelling shortwave radiation, interpolated value
+             REAL(wp) ::  sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value   
 
              IF ( surf%ns < 1 )  RETURN
+!
+!--          level-of-detail = 1. Note, here it must be distinguished between
+!--          averaged radiation and non-averaged radiation for the upwelling
+!--          fluxes. 
+             IF ( rad_sw_in_f%lod == 1 )  THEN
+             
+                sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var1d(tm)            &
+                                   + fac_dt * rad_sw_in_f%var1d(t)
+                                         
+                lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var1d(tm)            &
+                                   + fac_dt * rad_lw_in_f%var1d(t)
+!
+!--             Limit shortwave incoming radiation to positive values, in order 
+!--             to overcome possible observation errors.
+                sw_in = MAX( 0.0_wp, sw_in )
+                sw_in = MERGE( sw_in, 0.0_wp, sun_up )
                           
-             surf%rad_sw_in = sw_in                                         
-             surf%rad_lw_in = lw_in
-
-             IF ( average_radiation )  THEN
+                surf%rad_sw_in = sw_in                                         
+                surf%rad_lw_in = lw_in
              
-                surf%rad_sw_out = albedo_urb * surf%rad_sw_in
+                IF ( average_radiation )  THEN
+                   surf%rad_sw_out = albedo_urb * surf%rad_sw_in
+                   
+                   surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4  &
+                                  + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in
+                   
+                   surf%rad_net = surf%rad_sw_in - surf%rad_sw_out             &
+                                + surf%rad_lw_in - surf%rad_lw_out
+                   
+                   surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb          &
+                                                     * sigma_sb                &
+                                                     * t_rad_urb**3
+                ELSE
+                   DO  m = 1, surf%ns
+                      k = surf%k(m)
+                      surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m)  *      &
+                                             surf%albedo(ind_veg_wall,m)       &
+                                           + surf%frac(ind_pav_green,m) *      &
+                                             surf%albedo(ind_pav_green,m)      &
+                                           + surf%frac(ind_wat_win,m)   *      &
+                                             surf%albedo(ind_wat_win,m) )      &
+                                           * surf%rad_sw_in(m)
+                   
+                      surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m)  *      &
+                                             surf%emissivity(ind_veg_wall,m)   &
+                                           + surf%frac(ind_pav_green,m) *      &
+                                             surf%emissivity(ind_pav_green,m)  &
+                                           + surf%frac(ind_wat_win,m)   *      &
+                                             surf%emissivity(ind_wat_win,m)    &
+                                           )                                   &
+                                           * sigma_sb                          &
+                                           * ( surf%pt_surface(m) * exner(k) )**4
+                   
+                      surf%rad_lw_out_change_0(m) =                            &
+                                         ( surf%frac(ind_veg_wall,m)  *        &
+                                           surf%emissivity(ind_veg_wall,m)     &
+                                         + surf%frac(ind_pav_green,m) *        &
+                                           surf%emissivity(ind_pav_green,m)    &
+                                         + surf%frac(ind_wat_win,m)   *        &
+                                           surf%emissivity(ind_wat_win,m)      &
+                                         ) * 4.0_wp * sigma_sb                 &
+                                         * ( surf%pt_surface(m) * exner(k) )**3
+                   ENDDO
                 
-                surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4     &
-                                 + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in
-                
-                surf%rad_net = surf%rad_sw_in - surf%rad_sw_out                &
-                             + surf%rad_lw_in - surf%rad_lw_out
-                
-                surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb  &
-                                                  * t_rad_urb**3
-!
-!--          Calculate upwelling radiation fluxes and net radiation for each 
-!--          surface element depending on its properties.
+                ENDIF
+!
+!--             If diffuse shortwave radiation is available, store it on 
+!--             the respective files.
+                IF ( rad_sw_in_dif_f%from_file )  THEN 
+                   sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var1d(tm)  &
+                                         + fac_dt * rad_sw_in_dif_f%var1d(t)
+                                         
+                   IF ( ALLOCATED( rad_sw_in_diff ) )  rad_sw_in_diff = sw_in_dif
+                   IF ( ALLOCATED( rad_sw_in_dir  ) )  rad_sw_in_dir  = sw_in  &
+                                                                    - sw_in_dif
+!              
+!--                Diffuse longwave radiation equals the total downwelling 
+!--                longwave radiation
+                   IF ( ALLOCATED( rad_lw_in_diff ) )  rad_lw_in_diff = lw_in
+                ENDIF
+!
+!--          level-of-detail = 2
              ELSE
 
@@ -2801,4 +2957,17 @@
                    j = surf%j(m)
                    k = surf%k(m)
+                   
+                   surf%rad_sw_in(m) = ( 1.0_wp - fac_dt )                     &
+                                            * rad_sw_in_f%var3d(tm,j,i)        &
+                                   + fac_dt * rad_sw_in_f%var3d(t,j,i)                                   
+!
+!--                Limit shortwave incoming radiation to positive values, in 
+!--                order to overcome possible observation errors.
+                   surf%rad_sw_in(m) = MAX( 0.0_wp, surf%rad_sw_in(m) )
+                   surf%rad_sw_in(m) = MERGE( surf%rad_sw_in(m), 0.0_wp, sun_up )
+                                         
+                   surf%rad_lw_in(m) = ( 1.0_wp - fac_dt )                     &
+                                            * rad_lw_in_f%var3d(tm,j,i)        &
+                                   + fac_dt * rad_lw_in_f%var3d(t,j,i) 
 !
 !--                Weighted average according to surface fraction. 
@@ -2833,20 +3002,39 @@
                    surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)    &
                                    + surf%rad_lw_in(m) - surf%rad_lw_out(m)
+!
+!--                If diffuse shortwave radiation is available, store it on 
+!--                the respective files. 
+                   IF ( rad_sw_in_dif_f%from_file )  THEN
+                      IF ( ALLOCATED( rad_sw_in_diff ) )                       &
+                         rad_sw_in_diff(j,i) = ( 1.0_wp - fac_dt )             &
+                                              * rad_sw_in_dif_f%var3d(tm,j,i)  &
+                                     + fac_dt * rad_sw_in_dif_f%var3d(t,j,i)
+!
+!--                   dir = sw_in - sw_in_dif. 
+                      IF ( ALLOCATED( rad_sw_in_dir  ) )                       &
+                         rad_sw_in_dir(j,i)  = surf%rad_sw_in(m) -             &
+                                               rad_sw_in_diff(j,i)
+!                 
+!--                   Diffuse longwave radiation equals the total downwelling 
+!--                   longwave radiation
+                      IF ( ALLOCATED( rad_lw_in_diff ) )                       &
+                         rad_lw_in_diff = surf%rad_lw_in(m)
+                   ENDIF
 
                 ENDDO
 
              ENDIF
-                
-!
-!--          Fill out values in radiation arrays
+!
+!--          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
-             
-             rad_sw_in(0,:,:) = sw_in
-             rad_lw_in(0,:,:) = lw_in
  
           END SUBROUTINE radiation_external_surf
